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ABSTRACT 


Fog  situations  are  a hazard  to  navigation  whether  they  involve  highways,  har- 
bors, airports,  or  sea  lanes.  Prediction  and  modification  are  therefore  necessary 
for  the  safety  of  travel  and  for  operational  economy.  Field  studies  conducted  for 
this  purpose  indicate  that  many  different  fog  situations  result  from  a variety  of 
natural  processes,  which  must  be  systematically  investigated  so  that  the  basic 
causes  for  their  development  can  be  understood  and  possibly  controlled. 

The  systems  study  of  marine  fog  situations  that  has  been  conducted  has  taken 
into  consideration  the  dynamical  interaction  between  the  synoptical  variations  of 
meteorology  and  the  microphysical  phenomena  in  cloud  physics.  The  resulting 
numerical  program  simulates  the  transport  of  heat,  mass,  and  momentum  in  the 
atmospheric  boundary  layer  where  condensation/evaporation,  thermal  radiation, 
collision  coalescence,  and  turbulent  advection  take  place. 

Although  this  program  is  fairly  complete,  experimental  data"  of  the  moisture 
supply  process  affected  by  turbulent  winds  at  the  air-sea  interface  are  needed  to 
establish  a relaistic  model  for  systems  analysis.  Because  of  this  deficiency  in  the 
current  model,  the  radiation  estimate,  which  is  closely  related  to  the  water  vapor 
and  liquid  water  concentrations  at  various  heights  in  the  atmospheric  boundary  layer, 
can  only  be  considered  as  "qualitative".  When  available  field  data  are  compared  to 
predicted  fog  situations,  it  is  found  that  the  boundary  layer  model  can  qualitatively 
simulate  fog  situations  over  an  open  ocean,  but  the  accurate  prediction  of  a fog 
event  at  a specific  time  and  place  requires  a working  knowledge  of  the  moisture 
supply  process  occasioned  by  air-wave  interactions.  Until  this  knowledge  is  acquired, 
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the  "qualitative"  models  currently  being  used  in  systems  analysis  will  continue  to 
be  a source  of  inaccuracy. 

Fog  situations  are  known  to  occur  on  the  lee  side  of  some  topographical  ob- 
structions. The  wind  patterns  associated  with  this  type  of  fog  generally  have  dif- 
ferent directions  with  respect  to  the  prevailing  wind.  A hypothesis  is  being  formulated 
to  explain  such  situations  on  the  basis  of  recirculating  flows.  Some  "strange"  fog 
conditions  have  been  logically  explained  by  this  hypothesis.  However,  a realistic 
understanding  of  all  marine  fog  situations  requires  the  continuous  effort  of  dynamic 
modeling  to  integrate  laboratory  experiments  and  field  measurements. 


1.  INTRODUCTION 


Fog  situations  are  a serious  hazard  to  everyday  life  because  extensive  navi- 
gation in  the  air,  on  land,  and  at  sea  is  an  integral  part  of  travel,  commerce,  and 
defense.  Thousands  of  lives  and  millions  of  dollars  are  lost  each  year  because  of 
man's  inability  to  predict  and  to  modify  undesirable  fog  situations. 

Experimental  modification  of  fog  situations  has  been  attempted  at  harbors 
and  airports  by  Walker  and  Fox  (1942-1946),  Silverman  and  Kunkel  (1970),  Magono 
(1972),  and  Weinstein  (1973),  but  their  limited  success  makes  it  difficult  to  justify 
the  enormous  cost  for  practical  operation.  It  is  evident  that  an  understanding  of 
the  formation  mechanism  of  various  fog  situations  is  necessary  if  a realistic  method 
is  to  be  developed  for  successful  prediction  and  effective  modification. 

2.  FOGS  AND  CLOUDS 

The  similarity  of  fogs  and  clouds  prompted  the  author  to  concentrate  his  early 
efforts  on  the  study  of  basic  cloud  physics  problems,  such  as  condensation  nuclei, 
super  saturation,  and  collision  coalescence.  Field  data  by  Leipper  (1968),  Mack 
etaj.  (1972,  1975),  and  Larson  (1976)  indicate  that  the  marine  environment  has 
numerous  condensation  nuclei  that  do  not  allow  supersaturation  to  take  place  at  the 
air-sea  interface.  In  other  words,  fog  will  form  as  soon  as  excessive  moisture  is 
available  at  the  surface  layer  of  the  marine  environment.  Most  likely,  the  sea 
salt  concentration  in  marine  aerosols  may  cause  fog  to  form  even  in  subsaturated 
conditions,  because  sea  salt  nuclei  are  known  to  grow  in  size  by  absorbing  moisture 
in  a subsaturated  environment  (Mason,  1971). 

The  phenomenon  of  collision  coalescence  of  cloud  drops  in  a rain  forming  pro- 
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cess  appears  to  be  responsible  for  fog  dissipation  by  hygroscopic  particle  seeding 
as  discussed  by  Tag  (1971).  A hydrodynamic  analysis  was  developed  by  Lin  and  Lee 
(1973,  1975,  1976a,  1976b)  in  their  study  of  the  collision  efficiency  of  free  falling 
drops  in  the  atmosphere.  In  this  analysis,  it  became  evident  that  wake  capture  is 
a very  important  mechanism,  which  increases  the  collision  efficiency  of  a collecting 
drop  and  a collector  drop  when  their  radius  ratios  are  greater  than  90  percent.  How- 
ever, collision  efficiency  becomes  significantly  less  effective  than  geometric  efficiency 
when  the  radius  ratios  are  less  than  30  percent.  Considerable  discussion  was 
generated  by  this  analysis,  because  earlier  published  statements  claimed  that  wake 
capture  was  not  reinvent  to  an  increase  in  collision  efficiency.  Appendix  A contains 
the  four  papers  in  which  the  phenomenon  is  analyzed.  The  authors'  findings  indicate 
that  it  is  possible  to  produce  rainfall  through  the  process  of  collision  coalescence 
among  growing  drops  when  the  condensation  growth  of  the  seeded  particles  in  a 
supersaturated  cloud  reaches  the  critical  size  of  approximately  20  radii.  How- 
ever, the  application  of  cloud  seeding  is  not  effective  in  marine  fog  dissipation, 
because  there  are  enough  nuclei  already  present  in  the  air-sea  interface. 

The  necessary  ingredients  are  condensation  nuclei  and  saturation  conditions  for 
both  fog  and  cloud  formation;  however,  the  mechanisms  for  supplying  these  ingre- 
dients are  quite  different.  Clouds  can  generally  be  formed  by  the  cooling  of  the 
available  moisture  at  high  altitudes,  a condition  which  can  be  simulated  in  a cloud 
chamber.  Fogs,  on  the  other  hand,  can  be  formed  by  the  additional  mechanism  of 
providing  moisture  from  a nearby  water  source.  This  condition  cannot  be  simulated 
in  a conventional  cloud  chamber.  Consequently,  it  is  necessary  to  use  a systems 
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study  to  understand  the  dynamic  interaction  that  takes  place  between  synoptical 
meteorology  and  microphysics.  Moreover,  the  formation  of  clouds  in  severe 
storms — one  of  the  few  remaining  catastrophic  agents  against  which  mankind  has 
no  practical  defense — can  only  be  apprehended  through  a systematically  developed 
study  of  the  dynamical  and  microphysical  interactions  in  the  atmosphere. 

3.  A SYSTEMS  STUDY  OF  MARINE  FOG  SITUATIONS 
During  an  examination  of  the  field  data  of  Taylor  (1917),  Stone  (1936),  George 
(1940),  Emmons  (1947),  and  Magono  (1972),  it  became  clear  that  most  fog  situations 
are  caused  by  the  interaction  between  the  dynamical  process  of  heat,  mass  and 
momentum  transfer,  and  the  microphysical  process  of  condensation/evaporation, 
turbulent  advection,  and  thermal  radiation.  The  derived  systems  study  of  the  marine 
fog  situation  is  illustrated  diagramatically  by  the  flow  chart  in  Figure  1.  The  input 
information  consists  of  the  boundary  condition  of  the  potential  fog  region  and  the 
initial  condition  of  the  meteorological  data,  which  can  be  obtained  from  nearby 
weather  stations  and/or  available  weather  satellites.  The  dynamic  equations,  which 
are  given  in  Figure  2,  are  used  to  analyze  the  synoptical  meteorological  conditions 
that  govern  the  temperature,  humidity,  and  wind  velocity  of  the  potential  fog  region. 
The  microphysical  phenomena,  which  are  outlined  in  Figure  3,  are  concerned  with 
the  interdependent  processes  of  condensation/vaporation,  turbulent  advection,  collision 
coalescence,  and  thermal  radiation.  Laboratory  experiments,  which  are  often  con- 
ducted in  a well  controlled  environment,  can  only  provide  the  empirical  relation  of  a 
microphysical  process  as  a function  of  other  meteorological  parameters  in  a non- 
interacting environment.  To  consider  any  dynamic  interaction  between  the  s>  aoptical 
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meteorological  conditions  and  some  related  microphysical  phenomena,  it  is  necessary 
to  include  all  microphysical  processes?  as  the  appropriate  driving  functions  of  the 
dynamic  equations.  The  output  gives  the  primary  results  of  wind  velocity,  air 
temperature,  dew  point,  and  heat  budget  as  well  as  the  liquid  water  and  water  vapor 
contents.  The  drived  information,  such  as  visibility  and  refraction  index,  can  be 
calculated  from  the  primary  results  at  a given  condition  of  concentration  and  com- 
position of  the  marine  aerosols.  Field  data  can  then  be  used  to  verify  each  analyti- 
cal procedure  to  ensure  the  validity  and  reliability  of  the  developed  program  for 
marine  fog  simulations.  When  the  program  is  in  a ready  state  for  practical  applica- 
tion, weather  station  data  can  then  be  used  for  forecasting  fog  situations  in  a given 
region.  If  the  predicted  fog  is  expected  to  be  hazardous  to  navigation,  possible 
modification  schemes  can  be  simulated  to  determine  their  effectiveness  both  scienti- 
fically and  economically.  The  success  of  such  a systems  study  for  marine  fog 
simulation  is  based  on  two  factors: 

1.  The  Reliability  and  Adequacy  of  Field  Data 

Field  data  have  to  be  accurately  recorded  at  stratigic  locations  during  the  critical 
period.  For  example,  at  each  fog  event,  it  is  necessary  to  collect  all  fog  data  both 
inside  and  outside  the  fog  bank  in  line  with  the  direction  of  the  wind.  The  formation 
process  can  be  studied  by  using  the  data  taken  from  the  leading  edge  (or  upwind  region) 
of  the  fog  bank.  The  dissipation  process  can  be  studied  by  using  the  trailing  edge  (or 
downwind  region)  data.  Field  data  taken  in  the  middle  of  a mature  fog  bank  are  valu- 
able for  studying  its  chemical  composition  but  cannot  provide  the  required  information 
for  an  understanding  of  the  dynamic  interactions  that  occur  during  formation  or  dissipa- 


tion.  On  the  other  hand,  field  data  taken  across  fog  boundaries  in  a direction  normal 
to  the  wind  consist  of  too  many  uncertainties  and  cannot  provide  the  primary  infor- 
mation that  is  needed  for  an  understanding  of  the  formation  and  dissipation  processes. 
2.  The  Ability  of  Obtaining  Realistic  Mathematical  Solutions 

Numerical  solutions  of  nonlinear,  simultaneous,  partial  differential  equations 
can  be  obtained  by  using  high  speed  digital  computers,  but  marine  fog  problems,  which 
are  compunded  by  turbulent  winds,  require  an  appropriate  closure  scheme  to  describe 
the  physical  problem  realistically.  Because  of  the  physical  and  mathematical  com- 
plexities, marine  fog  analyses  are  primarily  made  of  open  ocean  situations. 

4.  OPEN  OCEAN  FOG  SITUATIONS 
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Fog  situations  that  originate  over  the  ocean  and  are  advected  to  coastal  regions 
have  been  reported  in  England  by  Walker  and  Fox  (1942-1946),  in  the  United  States 
by  Leipper  (1968),  in  Japan  by  Magono  (1972),  and  in  many  other  coastal  regions  of 
the  world.  A 48-hour  fog  event,  which  occurred  along  the  coast  of  Massachusetts 
during  the  26-28  January  1975,  demonstrates  the  feasibility  of  predicting  fog  events 
on  the  basis  of  weather  station  information.  A number  of  weather  stations  are  loca- 
ted in  the  region,  as  shown  in  Figure  4.  The  wind  direction,  wind  speed  (in  knots), 
air  temperature  (in  °C),  dew  point  (in  °C),  and  visibility  (in  km)  are  shown  in  the 
subsequent  figures  for  all  the  reporting  stations.  The  average  air  temperature  in 
the  region  before  the  event  was  about  3°C.  Warm  wind  moved  into  the  region  from 


1 the  east  and  southeast  direction  at  approximately  10  knots  during  the  early  morning 


hours  of  26  January  1975  (Fig.  4a).  Heavy  fog  began  to  appear  at  the  1500  hour  (Fig. 
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point  temperature  was  lower  than  the  local  air  temperature.  As  long  as  the  warm 
humid  air  was  supplied  from  the  south  or  southeast  direction,  fog  persisted  in  spite 
of  substantial  dew  point  depressions  (Fig.  4c).  The  fog  began  to  dissipate  in  the 
inland  area,  when  cold  dry  air  moved  into  the  region  from  the  west  at  the  1200  hour 
of  28  January  1975  (Fig.  4d).  Fog  was  cleared  from  the  entire  region  at  the  1800 
hour  (Fig.  4e),  when  the  average  air  temperature  was  about  2°C  and  the  wind  moved 
at  15  knots  from  the  west  direction.  During  the  45  hours  of  the  heavy  fog  period, 
two  observations  need  to  be  emphasized:  1)  The  diurnal  effect  of  radiation  heating 
did  not  appear  to  play  any  significant  role  in  improving  visibility,  and  2)  significant 
dew  point  depressions  were  reported  by  all  the  weather  stations  where  the  dew 
point  data  were  recorded. 

Advection  fog  of  this  type  is  a major  hazard  to  navigation,  and  this  particular 
fog  event  makes  it  very  clear  that  synoptical  meteorology  has  to  be  an  integral 
part  of  the  investigation.  Mathematical  formulation  of  fog  analysis  was  suggested 
by  Rodhe  (1962),  Fisher  and  Caplan  (1963),  and  Zakharova  (1972).  The  turbulent 
nature  of  the  transport  process  of  the  atmospheric  air,  as  discussed  by  Lumley 
and  Panofsky  (1964)  and  Monin  and  Yaglom  (1971),  limited  the  author's  choice  of 
realistic  closure  schemes  for  analyzing  actual  fog  situations.  Lee  et  al.  (1974) 
developed  a boundary  layer  model  to  study  the  turbulent  transport  of  heat,  mass, 
and  momentum.  Pepper  and  Lee  (1975)  extended  this  method  to  investigate  advection 
fog  over  the  ocean.  The  air-sea  interface  can  be  treated  as  an  atmospheric  boundary 
layer  that  has  various  roughnesses  at  the  ocean  surface.  The  amount  of  moisture 
supply  from  the  ocean  depends  on  the  air-water  temperature  difference  as  well  as 
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on  the  spray  process  resulting  from  air-wave  interactions.  In  general,  fog  may 
form  dui’ing  one  of  the  following  situations: 

(1)  Warm  humid  air  can  be  cooled  to  a saturation  condition  while  it  is  being 
transported  by  turbulent  winds  to  a cold  region.  This  was  the  condition  under 
which  the  above-mentioned  Massachusetts  fog  was  formed  as  well  as  the  advection 
fog  observed  by  Magona  (1972)  at  the  south  shore  of  Hakkaido,  Japan.  Marine  fog 
of  this  type  constitutes  the  most  undesirable  hazard  to  navigation,  because  it  may 
last  several  days  if  the  synoptical  conditions  persist.  Magono  (1972)  was  able  to 
improve  visibility  from  150  m to  250  m for  five  minutes  by  burning  five  tons  of 
propane  along  half  of  the  runway  length  at  the  Chitose  Airport.  This  example  makes 
it  clear  why  there  is  a need  to  study  the  scientific  and  economic  feasibility  of  fog 
modification. 

(2)  Moisture  may  be  supplied  from  a warm  ocean  to  cooler  air  through  evapora- 
tion. Fog  is  formed  when  radiation  cooling  drops  the  air  temperature  to  a conden- 
sation condition.  Wind  accelerates  the  evaporation  process.  Turbulent  winds  often 
add  more  moisture  to  the  cool  air  by  evaporating  the  spray.  The  more  moisture 
there  is  in  the  air,  the  less  degree  of  radiation  cooling  is  needed  to  reach  saturation 
conditions.  This  type  of  fog  often  occurs  as  a nocturnal  phenomenon  as  observed  by 
George  (1940)  along  the  Gulf  of  Mexico,  by  Mack  et  al.  (1972)  along  the  California 
Coast,  and  by  Goodman  (1976)  in  San  Francisco,  California.  Radiation  cooling  is  a 
very  important  factor  in  the  formation  of  dense  fog  in  this  situation.  In  other  words, 
when  such  a moisture  supply  process  occurs  in  the  morning,  as  observed  by  Kikuehi 
(1964)  at  the  Ishikari  Bay  located  on  the  north  shore  of  Hokkaido,  Japan,  a thin  layer 
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of  the  condensed  moisture  may  soon  be  dissipated  by  the  radiation  heat  created 
by  the  thermal  energy  of  the  rising  sun. 

The  dynamics  involved  in  the  systems  study  of  predicting  fog  conditions  at  sea 
are  relatively  well  known,  but  the  microphysical  phenomenon  of  moisture  supply 
at  the  air-sea  interface  is  not  completely  understood,  especially  when  turbulent 
winds  are  involved  in  the  process.  Moreover,  the  effect  of  radiation  cooling  can 
be  correctly  evaluated  only  when  the  amounts  of  water  vapor  and/or  liquid  water 
in  the  radiating  layer  are  accurately  known.  The  adjacent  water  is  the  source  for 
the  moisture  that  creates  a marine  fog  situation.  This  is  quite  different  from  the 
situation  in  which  a cloud  is  formed  at  a relatively  high  altitude.  Laboratory 
experiments  cannot  ignore  the  problem  of  moisture  supply,  simply  because  the 
condition  is  unfamiliar  to  cloud  physicists.  Field  data  have  to  be  collected  to  pro- 
vide the  maximum  amount  of  information  that  can  be  attained  in  the  direction  of 
the  wind  from  both  outside  and  inside  of  a fog  bank. 

Detailed  discussions  of  the  boundary  layer  method  that  is  used  for  predicting 
marine  fog  situations  at  sea  are  contained  in  the  technical  papers  included  in 
Appendix  B.  A computer  program  for  the  method  is  available  upon  request.  Be- 
cause the  evaluation  of  a moisture  supply  cannot  be  made  with  certainty,  the  developed 
program  can  be  used  only  for  a qualitative  analysis  at  the  present.  The  data  that  are 
currently  available  from  laboratory  experiments  and  field  observations  should  be 
used  to  improve  the  microphysical  models  so  tha^  the  current  program  can  be 
developed  into  a practical  tool  for  fog  prediction.  Moreover,  it  must  be  noted  that 
the  boundary  layer  method  cannot  be  applied  to  fog  situations  in  which  the  topography 
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and/or  weather  fronts  cause  recirculating  flows. 


5.  FOG  FORMATION  IN  RECIRCULATING  FLOW  REGIONS 
Fog  may  form  along  shorelines  or  in  harbors  where  the  wind  recirculates  as 
a result  of  the  topography.  The  importance  of  recirculating  flow  as  a natural  pro- 
cess was  first  suggested  by  von  Karman  (1934).  In  his  autobiography  (von  Karman, 
1967),  he  stated  that  the  concept  of  the  Karman  vortex  could  not  be  accepted  by  the 
fluid  mechanics  community  until  the  collapse  of  the  Tacoma-Narrow  bridge  could  not 
be  explained  by  any  other  known  phenomenon.  Consequently,  the  reluctance  of  the 
meteorological  community  (especially  in  the  restricted  area  of  cloud  physics)  to 
investigate  the  hypothesis  that  fog  may  form  in  the  region  of  recirculating  flows  is 
understandable.  However,  it  would  be  irresponsible  to  discard  the  facts,  because 
they  cannot  be  satisfactorily  explained  by  the  current  knowledge  of  cloud  physicists. 
The  facts  are  as  follows: 

(1)  The  fog  situation  observed  by  Mack  et  al.  (1975)  during  the  expedition  of  the 
R/V  Acania  in  the  region  downwind  of  Cape  Mendocino,  California,  on  25-27 
August  1974  was  considered  "strange”  (Fig.  5). 

(2)  The  fog  situation  observed  by  Gathman  (1976)  and  Larson  (1976)  during  the 
expedition  of  the  USNS  Hays  in  the  region  downwind  of  Nova  Scotia,  Canada, 
on  2-3  August  1975  was  never  explained  (Fig.  6 and  Table  1). 

These  two  fog  situations,  which  differ  from  most  other  reported  fogs,  are  similar 
in  the  following  respects:  1)  In  both,  the  prevailing  winds  were  relatively  strong. 

2)  The  wind  directions  in  the  fog  regions  were  opposite  to  those  of  the  prevailing 
winds.  3)  Both  fog  regions  are  located  on  the  lee  side  of  topographic  obstructions. 


I 
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A complete  analysis  of  the  flow  characteristics  in  a recirculating  flow  region 


requires  the  solution  of  the  complete  set  of  equations  shown  in  Figure  2 without  the 
use  of  the  boundary  layer  assumption  outlined  in  the  technical  papers  contained  in 
Appendix  B.  There  are  two  problems.  One  is  to  develop  a numerical  method  for 
stable  solutions  of  a system  of  nonlinear  elliptical  differential  equations,  and  the 
other  is  to  establish  a realistic  closure  scheme  in  which  all  the  tensor  components 
of  turbulence  correlation  for  velocities,  temperatures,  and  water  vapor  and  liquid 
water  concentrations  can  be  considered.  The  author’s  effort  to  resolve  the  mathe- 
matical problem  was  discussed  by  Lin  et  a_l.  (1976).  This  paper  is  included  in 
Appendix  C.  The  Strongly  Implicit  Procedure  (SIP)  appears  to  be  the  most  efficient 
method  for  analyzing  recirculating  flow  when  the  turbulence  closure  scheme  is 
replaced  by  an  equivalent  laminar  model.  A concept  of  a turbulent  recirculating 
flow  can  be  developed  by  using  the  SIP  method  of  Appendix  C and  the  turbulence 
closure  scheme  of  Appendix  B which  still  needs  to  be  perfected  for  use  in  recirculating 
flow  regions.  Ability  to  predict  and/or  modify  fog  situations  in  harbors  requires  a 
basic  understanding  of  the  recirculating  flow  affected  by  topography.  This  knowledge 
is  not  currently  available. 

6.  CONCLUSIONS 


*•  1.  Fog  can  be  formed  by  either  decreasing  the  air  temperature  or  increasing  the 

air  humidity  to  a saturation  condition.  Once  the  saturation  condition  is  reached, 


that  have  not  been  realistically  simulated  in  conventional  cloud  chambers. 

3.  Prediction  and  modification  of  marine  fog  situations  cannot  be  accomplished 
by  considering  condensation  from  the  microphysical  viewpoint  alone,  because 
the  interaction  between  dynamics  and  microphysics  requires  the  synoptical 
information  of  dynamical  meteorology  and  oceanography  as  well  as  the  micro- 
physical phenomena  of  cloud  physics. 

4.  Because  there  are  many  possible  combinations  of  natural  phenomena  in  fog 
formations,  a systems  study  is  necessary  to  obtain  a general  perspective  in 
analyzing  the  fog  data  that  are  gathered  from  various  locations  at  different 
times. 

5.  A summary  of  ten  reported  fog  situations  is  given  in  Table  2.  Eight  cases  can 
be  simulated  by  using  the  boundary  layer  model.  The  other  two  cases,  how- 
ever, can  only  be  explained  by  using  the  hypothesis  of  recirculating  flow. 
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Def.  turbulence  kinetic  energy , 


MICROPHYSICAL  PROCESSES 
Eddy  Transport  ( Pepper  ft  Lee  , 1975) 

Km  = fJL-pvSu' / =0.3^Q/l-gj-| 

Kh  = K -yC>wV  / - Km  /crh 

Kc  = D —pVJC'  / ^ - Km/  Cr 
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Kq  = -/)W(^-+  Q )/f|"  = Km/crq 

Radiation  Effect  ( Mack, et  al,  1972 *,  Barker, 1973) 
^ - \.6PcrT^Kk(fEXp[-\.eK)J2'dIL  dz] 


bz 


F#Zr)  = -(Tlf[  I — €(Zr,Zt  )j-  dZ 

rU)=  °"C[l  -€(Zr,0)>/JcrT4^^dZ 
Evaporation  Process  (Dankerwerts,  1951  ;Knudson&  Katz,  1958) 
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Average  contact  time  t=  k, 
Condensation  Rate  (McDonald,  1963;  Asai,l965) 

C. 

Precipitation  Velocity  (Mack,et  al,l972) 

, m 2/3 

Vd  = 5.3  X 10s  ( w / N ) 


Figure  3 


Microphysical  Equations 
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Fig.  4a,  Massachusetts  Fog  0100  HR.  26  January  1976 


Air  Temp.  6°C  Visibility  0 km 


Fig.  4b,  Massachusetts  Fog  1500  HR.  26  January  1Q76 
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Fig.  4d,  Massachusetts  Fog  1200  HR.  29  January  1976 
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f Abstract  Transient  state  solutions  of  the  Navier  Stokes  equations  were  obtained  for  incom- 

pressible flow  around  a sphere  accelerating  from  zero  initial  velocity  to  its  terminal  free  falling 
velocity  By  assuming  rotational  symmetry  about  the  axis  in  the  direction  of  motion,  the  Navier- 
Stokes  equations  and  the  continuity  equation  were  simplified  in  terms  of  vorticity  and  stream 
function.  The  instantaneous  acceleration  of  the  falling  sphere  was  calculated  by  considering 
the  difference  between  the  gravitational  force  and  the  drag  force  in  a transient  state  A set  of 
implicit  finite  difference  equations  was  developed.  In  order  to  obtain  accurate  information 
around  the  body,  an  exponential  transformation  along  the  radial  direction  was  used  to  provide 
liner  meshes  in  the  vicinity  of  the  surface  of  the  sphere.  The  vorticity  equation  was  solved  by 
an  alternating  direction  implicit  (ADI)  method  while  the  stream  function  equation  was  solved 
by  a successive  over-relaxation  (SOR)  method.  Simultaneous  solutions  were  obtained  Transient 
state  solutions  were  compared  with  steady  state  solutions  for  Reynolds  numbers  up  to  300 
Separations  first  occurred  at  a Reynolds  number  20  for  steady  state  flows  and  at  Reynolds 
numbers  22  4h  and  28-24  for  transient  state  flows  with  terminal  Reynolds  numbers  of  100  and 
300,  respectively.  Separation  angles,  sizes  of  separation  regions,  and  drag  coefficients  were 
calculated  for  both  steady  and  unsteady  states.  Good  agreement  was  obtained  with  existing 
experimental  data  in  the  steady  state. 


NOMENCLATURE 
a lattice  spacing  in  radial  direction 
A location  of  grid  point  in  the  positive  z direction 
h lattice  spacing  in  angular  direction 
B location  of  grid  point  in  the  positive  8 direction 
C location  of  grid  point  in  the  negative  z direction 
C„  total  drag  coefficient 
C'Dr  drag  coefficient  due  to  skin  friction 
C,lr  drag  coefficient  due  to  surface  pressure 
</  diameter  of  sphere 

/)  location  of  grid  point  in  the  negative  8 direction 
g gravitational  acceleration 
I wake  length 

0 location  of  grid  point,  the  origin 

p static  pressure  at  surface  of  sphere 
/’  static  pressure  in  free  stream 
r radial  distance 

r , radial  distance  of  outer  boundary 
R radius  of  sphere 
(Re),  local  Reynolds  number 
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/'  R terminal  Reynolds  number 
V terminal  free  stream  velocity 
V 5 local  free  stream  velocity 
; dimensionless  radial  coordinate 
0 angular  spherical  coordinate 
H*  stream  function 
10  vorticity 
/i  dynamic  viscosity 

[ i kinematic  viscosity 

p density  of  the  medium 
p,  density  of  the  sphere 

INTRODUCTION 

Current  interest  in  the  motions  of  small  particles  in  the  atmosphere  has  lead  to  a large 
number  of  extensive  investigations  of  flow  around  spherical  bodies.  Approximate  solutions 
[ for  steady  state  laminar  flow  at  very  low  Reynolds  numbers  (Re  < 1)  have  been  obtained 

by  Stokes[l]  and  Oseen[2].  Stokes  assumed  that  the  inertia  force  was  negligible.  Oseen  used 
a small  perturbation  technique  to  linearize  the  equation  of  motion.  The  exact  solution  to 
i Oseen's  linearized  equation  was  obtained  by  Goldstein[3].  Based  on  Goldstein’s  analysis, 

the  flow  pattern  around  a spherical  body  was  then  calculated  in  detail  by  Tomotika  and 
Aoi[4]  and  Pearcey  and  McHugh[5],  Improvement  of  these  analyses  was  made  by  Proudman 
and  Pearson[6]  for  the  purpose  of  covering  higher  Reynolds  number  regions.  However,  from 
the  experimental  work  of  Maxworthy[7],  this  improvement  was  found  to  be  valid  only  for 
Reynolds  numbers  below  1-3.  In  order  to  analyze  the  flow  patterns  in  the  steady  state, 
where  the  inertia  force  cannot  be  ignored,  numerical  solutions  of  the  nonlinear  Navier- 
| Stokes  equations  were  obtained  by  Jenson[8]  for  Reynolds  numbers  up  to  40.  With  modern 

I • computers,  Hamielec  el  al. [9]  improved  Jenson's  analysis  for  Reynolds  numbers  up  to  100 

and  Le  Clair  et  «/. [10]  refined  the  numerical  computation  to  cover  the  Reynolds  number 
I region  up  to  400.  Le  Clair  et  al.’ s results  agree  well  with  the  empirical  relations  of  Pruppacher 

and  Steinberger[l  1]  and  Beard  and  Pruppacher[12].  Numerical  solution  of  the  unsteady  state 
Navier  Stokes  equations  was  obtained  by  Rimon  and  Cheng[l3].  However,  when  the  steady 
f state  was  approached.  Rimon  and  Cheng's  results  differed  from  those  of  Le  Clair  el  al.  for 

separation  angle  in  the  Reynolds  number  range  between  10  and  20,  for  surface  pressure 
distributions  for  ail  of  the  cases  studied  l Re  < 300),  and  for  surface  vorticity  distributions 
[•  in  the  separated  flow  region  for  Reynolds  numbers  greater  than  100.  Rimon  and  Cheng's 

method  was  used  by  Shafrir  and  Tzvi[l4],  with  an  improved  boundary  condition,  for 
terminal  Reynolds  numbers  up  to  104  Shafrir  and  Tzvi's  solutions  gave  closer  agreement 
with  Le  Clair  et  a/.'s  results  than  those  of  Rimon  and  Cheng.  Experimental  data  ofTaneda[l  5] 
E>  *.  also  seems  to  support  Le  Clair  et  al.' s findings.  This  study  was  initiated  to  develop  a new 

method  for  transient  state  analysis  to  evaluate  the  flow  pattern  around  an  accelerating 
particle  of  spherical  shape. 

■ *■  * 

| THEORETICAL  ANALYSIS 

f „ , The  governing  equations  for  an  incompressible  flow  around  a spherical  body  with  rota- 

tional  symmetry  in  the  flow  direction  may  be  written  as  follows: 

= r sin  0 to  ( I ) 
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T.R.  terminal  Reynolds  number 
V terminal  tree  stream  velocity 
i\  local  free  stream  velocity 
; dimensionless  radial  coordinate 
0 angular  spherical  coordinate 
"P  stream  function 
w vorticity 
/i  dynamic  viscosity 
v kinematic  viscosity 
p density  of  the  medium 
p,  density  of  the  sphere 

INTRODUCTION 

Current  interest  in  the  motions  of  small  particles  in  the  atmosphere  has  lead  to  a large 
number  of  extensive  investigations  of  flow  around  spherical  bodies.  Approximate  solutions 
for  steady  state  laminar  flow  at  very  low  Reynolds  numbers  (Re  < I)  have  been  obtained 
by  Stokesfl]  and  Oseen[2).  Stokes  assumed  that  the  inertia  force  was  negligible.  Oseen  used 
a small  perturbation  technique  to  linearize  the  equation  of  motion.  The  exact  solution  to 
Oseen's  linearized  equation  was  obtained  by  Goldstein[3].  Based  on  Goldstein's  analysis, 
the  flow  pattern  around  a spherical  body  was  then  calculated  in  detail  by  Tomotika  and 
Aoi[4J  and  Pearcey  and  McHugh[5],  Improvement  of  these  analyses  was  made  by  Proudman 
and  Pearson  [6]  for  the  purpose  of  covering  higher  Reynolds  number  regions.  However,  from 
the  experimental  work  of  Maxworthy[7],  this  improvement  was  found  to  be  valid  only  for 
Reynolds  numbers  below  1-3.  In  order  to  analyze  the  flow  patterns  in  the  steady  state, 
where  the  inertia  force  cannot  be  ignored,  numerical  solutions  of  the  nonlinear  Navier  - 
Stokes  equations  were  obtained  by  Jenson[8]  for  Reynolds  numbers  up  to  40.  With  modern 
computers,  Hamielec  or  a/.[9]  improved  Jenson’s  analysis  for  Reynolds  numbers  up  to  100 
and  Le  Clair  el  u/.[IO]  refined  the  numerical  computation  to  cover  the  Reynolds  number 
region  up  to  400.  Le  Clair  el  al.'s  results  agree  well  with  the  empirical  relations  of  Pruppacher 
and  Steinberger[l  I ] and  Beard  and  Pruppacfier[l  2).  Numerical  solution  of  the  unsteady  state 
Navier  Stokes  equations  was  obtained  by  Rimon  and  C'heng[l  3],  However,  when  the  steady 
state  was  approached,  Rimon  and  Cheng's  results  differed  from  those  of  Le  Clair  el  al  for 
separation  angle  in  the  Reynolds  number  range  between  10  and  20,  for  surface  pressure 
distributions  for  ail  of  the  cases  studied  I Re  < 300),  and  for  surface  vorticity  distributions 
in  the  separated  flow  region  for  Reynolds  numbers  greater  than  100.  Rimon  and  Cheng's 
method  was  used  by  Shalrir  and  Tzvi[  1 4],  with  an  improved  boundary  condition,  for 
terminal  Reynolds  numbers  up  to  104  Shafrir  and  Tzvi's  solutions  gave  closer  agreement 
with  LeClair  el  ul.\ results  than  those  of  Rimon  and  Cheng,  experimental  data  of  Taneda[l  5] 
also  seems  to  support  Le  Clair  el  ill's  findings.  This  study  was  initiated  to  develop  a new 
method  for  transient  state  analysis  to  evaluate  the  flow  pattern  around  an  accelerating 
particle  of  spherical  shape. 

I H FOR  I TIE  At.  ANALYSIS 

The  governing  equations  for  an  incompressible  flow  around  a spherical  body  with  rota- 
tional symmetry  in  the  flow  direction  may  be  written  as  follows: 

r sin  0 i,i  (I) 
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where  F = • 


cu 


e'  sin  (I 
G = me'  sin  0 

r.i 


t2  = k->-k  + s'ne7o{™ok) 


Equations  (7),  (8)  and  (9)  may  be  solved  simultaneously  for  the  stream  function.  S',  the 
vorticity,  co,  and  the  instantaneous  Reynolds  number.  ( Re ),,  in  terms  of  the  independent 
variables  of  /,  r and  I). 

FINITE  DIFFERENCE  EQUATIONS 

Numerical  solutions  can  be  obtained  by  using  central  difference  and  forward  difference 
approximations  for  the  space  and  time  coordinates,  respectively.  The  finite  difference 
equations  may  be  written  as  follows: 

I . Stream  function 

The  stream  function  equation  (7)  may  be  written  in  the  following  finite  difference  form, 
which  can  be  solved  by  a successive  over-relaxation  (SOR)  method: 

+ HV8  + TV  Ac  + 4V  A„  t T„;.0  - (me3r  sin  0)o  = 0 

whereA^ig-l) 


(10) 


'■  fji  - ■-“'"•I 

'""sic' 


with  a - Ac  and  b A 0.  The  subscripts  .4.  B.  C.  O and  0 refer  to  the  locations  of  the  grid 

points  in  space  as  shown  in  Tig  I 

2.  Vorticity 

The  vorticity  equation  (8)  may  be  written  as  two  finite  difference  equations  which  can  be 
solved  by  an  alternating  directional  implicit  (ADI)  scheme: 

(I)  The  first  half  time-step.  (n  + J).  is  given  in  the  0 direction  as 


Co  + C = 0 


(ID 


where  C»  = 


(Re), 


sin  (K 


I6u/>e  ' Mn  On  am  wo 

_ iRe)'  w •*'  -V  + V --H*  -1  - Sin(?n 
(6ahei’sinODl  A ' A sin  0oe: 


cot  0o 
2b 


I cot  0o 

P + 
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G Y - 2G0"  + Gc"  CA"  - Gc" 

^0" 

d(/G),' 

a2  2 a 

(«e), 

dr 

2 2 
;°  = A 'l  + bI^i 

c - - ^7  - ['♦v*1  - v,r'  + v,-  - 4VH/V  - f <"] 

Af  I babe 

\ 0/  - 2G0"  + 6'c"  G/  ~ GY  m0"  I d(/tf  ),1" 

sind0e*'  a2  2a  + (Re),  l dr 

(2)  The  second  half  time-step,  (n  + 1 ),  is  given  in  the  r direction  as 

"Y*'sa  + <V"£o  + wc"  'Zc  + £ =0 

(Kei  , 111 

where  (, . = --r— -! [MV  - MV  + MY  - My]  - -y 5 

3 16aheJ  sin  0o  1 " b nt  e2'  “ u2  2u 

S ^ ^ ^ [UJ  "+l  111  3 + I , UJ  « UJnl  ^ ^ 1 

Sc  - .310  7T  lTS  ~ ~D  ■*'  “»  " T0  * 2i  * o — I + T~ 

I6«he  sin  0o  e2z  “[a2  2 a 

2 2 

£0  ~ t. ; + ,,2-22 


At  a e 

2io0"* 


i = - [IV* ' - *PC-  + TY  - My][fs"  + * - VM] 

1 GYM-2G0”+>  + GY+*  gv4*-g0"+M 

— r — tz r; cot  (L 

sin  fl0e}  b 2 2h 

, rcKfteU" 

(Re),  [ dt 

The  superscripts  n./r  + j and  n + J designate  the  present,  the  first-half  and  the  second-half 
time  steps,  respectively.  The  term  (d(Re),/dt)"  denotes  the  slope  of  the  instantaneous 
Reynolds  number  at  the  present  time  step  and  (Re),  is  the  average  value  of  (Re)"  and 
(Re),"*'. 
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Fig  I Schematic  diagram  for  boundary  conditions  and  nomenclatures  used  in  finite  difference 

equations. 
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3 Instantaneous  Reynolds  number 

The  instantaneous  Reynolds  number  tan  be  calculated  by  using  two  successive  numeri- 
cal techniques  as  given  by  Conte[l6], 

(1)  The  Reynolds  number  at  the  (n  + l)th  time  step  may  be  estimated  by  the  Adams 
Bashforth  prediction  formula: 

(Re),"*1  = (Re)"  + ^ (55/*-  59/"-'  + 37/"‘2  -9/-3]  (13) 

(2)  The  estimation  may  then  be  refined  by  the  Adams-Moulton  correction  formula: 

(Re),"*1  = (Re),"  + ^[9/-*'  + I9/--5/"-'  + f-2]  (14) 

where  f " denote  the  slope  of  the  instantaneous  Reynolds  number  at  the  nth  time  step 

BOUNDARY  CONDITIONS 

The  initial  condition  is  given  such  that  the  sphere  begins  to  fall  when  the  surrounding 
medium  is  not  yet  disturbed: 

t=0:  H*  = 0,  io  = 0,  (Re),=0,  everywhere  (15) 

The  boundary  conditions  are  such  that,  for  all  time  intervals,  the  axis  of  symmetry 
remains  undisturbed;  the  surface  velocity  is  zero;  and  the  free  stream  denotes  a uniform 
flow  velocity  for  the  sphere  The  dimensionless  form  of  the  boundary  conditions,  as  shown 
schematically  in  Fig  I,  may  be  written  as. 


0 < l < oo : 

V =0, 

to  = 0, 

at  8 = 0 for  all  z 

= 0, 

to  — 0, 

at  8 = n for  all  z 

S'  =0, 

1 d2^ 
sin  f)  dz1 

at  z = 0 for  all  0 

4*  = Je2'  sin2  8, 

to  - 0, 

at  z -*  oo  for  all  8 

Numerical  solutions  may  then  be  obtained  by  solving  the  stream  function  equation,  the 
vorticity  equation  end  the  Reynolds  number  equation  simultaneously  for  a given  size  sphere 
falling  from  zero  initial  velocity  to  its  terminal  velocity.  No  special  problems  were  encoun- 
tered in  applying  the  ADI  method  to  the  present  problem,  either  in  the  interior  or  at  the 
boundaries  However,  the  amount  of  computer  time  employed,  while  modest  for  low 
Reynolds  number  (five  to  ten  minutes  on  the  C'DC  7600)  increases  substantially  at  the 
larger  Reynolds  numbers. 

NUMERII  AL  SOLUTIONS 

Solutions  for  stream  function  and  vorticity  were  obtained  for  both  steady  and  unsteady 
states  Figure  2 shows  the  steady  state  solutions  for  Reynolds  numbers  of  5.  20,  40.  100,  200 
and  3(X)  It  can  b;  seen  that  the  recirculation  region  begins  to  develop  at  Reynolds  number 
20  and  increases  its  size  as  the  Reynolds  number  increases  Moreover,  the  vorticity  gradient 
increases  very  rapidly  with  Reynolds  numbers  at  the  surface  of  the  sphere.  Figure  3 shows 
the  transient  state  solutions  for  a falling  sphere  at  terminal  Reynolds  numbers  IT  R lot  100 
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(Streamlines 
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Fig.  2.  Distributions  of  streamlines  and  vortices  around  a sphere  for  steady  state  flow 


and  300.  with  local  Reynolds  numbers  (Re),  of  5,  40  and  90  It  is  noted  that,  at  the  same 
local  Reynolds  number,  the  recirculating  region  at  T.R.  = 100  is  larger  than  that  at 
T.R.  = 300.  Comparison  between  the  steady  slate  and  the  transient  state  solutions  is  shown 
in  Fig.  4 It  is  evident  that  a smaller  recirculating  region  occurs  in  accelerating  flow  in  com- 
parison with  the  fully  developed  steady  state  flow  at  the  same  local  Reynolds  number 


DR A(i  EVALUATION 

The  drag  force  of  a falling  sphere  consists  of  the  skin  friction  drag  and  pressure  drag 
The  skin  friction  drag  coefficient  can  be  calculated  through  the  vorticity  as: 
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Fig.  3.  Distributions  of  streamlines  and  vortices  around  a sphere  for  transient  state  floss 

The  pressure  drag  coefficient  can  be  obtained  through  the  vorticity  and  stream  function 
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Figures  5 and  6 show  the  vorticity  and  pressure  distributions,  respectively,  at  the  surface 
of  the  sphere  at  various  local  Reynolds  numbers  for  both  steady  state  and  unsteady  state 
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Fig  4 Comparisons  of  streamlines  and  vortices  around  a sphere  between  steady  state  flow  and 

transient  state  flow. 


flows.  The  total  drag  coefficient,  as  shown  in  Fig.  7,  was  calculated  by  summing  the  skin 
friction  drag  and  the  pressure  drag  Jt  can  be  seen  that  the  drag  coefficient  in  the  steady  state 
agrees  well  with  existing  experimental  data  as  given  by  Schlichting[  17].  An  empirical 
relation  for  the  steady  state  drag  coefficient  is  proposed  as: 


for  Reynolds  numbers  up  to  1000.  Available  formulae  for  calculating  the  drag  coefficient 
are  listed  in  Table  1 for  reference.  Transient  state  drag  coefficients  were  also  calculated  for 
terminal  Reynolds  numbers  of  100  and  300  and  are  shown  by  broken  lines  in  Fig.  7. 
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Fig  5 Vorticity  distributions  on  the  surface  of  a sphere. 


Transient  State 


Fig.  6.  Pressure  distributions  on  the  surface  of  a sphere. 


FLOW  SEPARATION 

Comparisons  with  experimental  data  for  separation  angle  and  wake  length  are  shown  in 
Figs.  8 and  9,  respectively.  The  steady  state  solutions  indicate  that  flow  separation  occurs 
at  a Reynolds  number  of  20.  The  unsteady  solutions  indicate  that  separation  occurs  at 
Reynolds  numbers  22-46  and  28  24  for  terminal  Reynolds  numbers  of  100  and  .100.  respec- 
tively. Taneda  reported  that  oscillations  of  the  separated  fiow  region  were  observed  for 
terminal  Reynolds  numbers  greater  than  130.  Preliminary  results  of  the  numerical  solutions 
also  indicated  some  oscillatory  behavior;  however  this  was  later  confirmed  to  be  the  result 
of  numerical  instability. 
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Fig.  7.  Drag  coefficient  vs.  Reynolds  number  for  a sphere. 


90* 


80* 


70* 


60* 


□)  50 


a 

O) 

CO 


30  *L 


20* 


10* 


Experiment 
Taneda(1 956) 
d = 1 9 82  mn 
d = 1 5 08  mm 
Steady  State  Analysis 
Proudman  & Pearson 
(1957) 

Jenson  (1959) 

Hamielec  (1  967) 

Le  Clair  ( 1 970) 
Unsteady  State  Analysis 


r 5^pn  &Chenjj(1969) 


Shafr.  &Tzvi  (1971 
. P'esent  Results 
40*|-  S'eady 
U isteady 


(T  R 100) 
(T  R 300) 


$ 


r_ 


Reynolds  number,  [Re J, 
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EFFECTS  OF  STEP  SIZE  AND  OUTt  K BOUNDARY  ON 
NUMERICAL  CALCULATIONS 

Numerical  instability  usually  leads  to  divergence  of  the  numerical  solution.  However, 
under  certain  conditions  an  unstable  or  oscillatory  solution  may  be  misinterpreted  as 
describing  a physical  phenomenon  The  transient  state  anal;  for  separated  flow  around 
a sphere  is  one  such  example.  In  this  case,  the  oscillation  is  largely  due  to  the  choice  of 
step  size 

I he  grid  sizes  Ac,  AO,  At  are  considered  to  be  sufficiently  small,  when  a decrease  of  the 
step  size  no  longer  affects  the  solution,  viz.,  to  some  prescribed  tolerance.  Table  2 shows  the 
mesh  sizes  chosen  in  this  study.  It  is  necessary  to  point  out  that  the  required  mesh  size  can 
be  a function  of  the  terminal  Reynolds  number  T he  mesh  sizes  of  Ar  = OT  and  AO  = 6 
were  originally  chosen  for  all  Reynolds  numbers  and  stable  solutions  were  reached  for  all 
cases  with  terminal  Reynolds  numbers  less  than  100.  However,  at  terminal  Reynolds 
numbers  of  200  and  300.  the  calculated  flow  patterns  in  the  wake  region  of  the  sphere 
fluctuated  within  certain  finite  limits.  The  dashed  lines  in  Fig.  10  depict  the  variation  of 
wake  length  and  separation  angle,  with  respect  to  time,  at  various  terminal  Reynolds 
numbers  ranging  between  40  and  300.  Experimental  data  obtained  by  Taneda[!7j  also 
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Table  2.  Step  sizes  and  outer  boundaries  used  in  numerical  computations 
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showed  such  an  oscillatory  motion  at  terminal  Reynolds  numbers  larger  than  130.  In  order 
to  have  a reasonable  assurance  that  this  phenomenon  was  not  caused  by  numerical  insta- 
bility, the  mesh  size  was  reduced  to  A z = 0 03  and  A 0 = 3°.  The  solid  lines  in  Fig.  10  show 
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Fig.  10.  Effects  of  step  sizes  on  numerical  instability  at  steady  state. 

that  the  solutions  are  stable  even  for  the  cases  with  terminal  Reynolds  numbers  of  200  and 
300  if  the  mesh  sizes  are  sufficiently  small.  It  is  therefore  necessary  to  conclude  that  the 
oscillatory  phenomenon  obtained  in  the  present  solution  with  larger  values  of  Ac  and  A 0 
is  due  to  numerical  instability. 
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The  question  of  whether  130  is  the  critical  Reynolds  number  in  a laminar  flow  around 
a sphere  cannot  be  answered  in  this  study;  the  analysis  assumes  rotational  symmetry  which 
may  present  the  development  of  a realistic  asymmetric  flow  pattern.  It  is  felt  that  the  critical 
Reynolds  number  can  only  be  determined  numerically  by  analyzing  the  complete  three 
dimensional  transient  flow  field  around  a spherical  body. 

The  accuracy  of  the  numerical  solution  is  also  affected  by  the  location  of  the  outer 
boundary  Theoretically,  the  limit  of  the  outer  boundary  is  at  infinity.  Unfortunately,  such 
a limit  cannot  be  used  in  numerical  computations.  Some  large  but  finite  distance  from  the 
sphere  therefore  has  to  be  assumed  The  distance  is  determined  in  such  a way  that  further 
increases  in  the  outer  boundary  do  not  affect  the  solution.  Table  3 indicates  the  effect  on  the 


Table  3.  Fffect  of  outer  boundary  on  steady  state  drag  calculations 
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solution  of  the  choice  of  the  outer  boundary  location.  As  the  Reynolds  number  increases,  the 
effect  of  the  outer  boundary,  on  the  accuracy  of  the  numerical  solution,  gradually  decreases 
while  the  effect  of  decreasing  step  size  becomes  more  important. 

CONCLUSION 

A transient  stale  analysis  has  been  formulated  for  separated  flow  around  a sphere.  The 
analytical  results  agree  well  with  available  experimental  data  and  other  analytic  solutions 
when  the  steady  state  is  approached.  The  transient  development  of  the  separation  region  is 
found  to  be  a function  of  the  local  Reynolds  number.  In  the  steady  state,  flow  separation 
occurs  at  a Reynolds  number  of  20.  In  the  transient,  flow  separation  occurs  at  local  Reynolds 
numbers  of  22  46  and  28  24  for  accelerating  spheres  with  terminal  Reynolds  numbers  of  100 
and  300.  respectively. 

The  drag  coefficient  of  a sphere  moving  at  a constant  velocity  is  usually  determined  from 
an  empirical  formula.  Based  on  the  comparison  between  analytical  and  experimental 
results,  a simple  relation  is  obtained  for  calculating  the  drag  coefficient  of  a sphere  moving 
uniformly  This  expression  applies  in  the  Reynolds  number  range  between  one  and  one 
thousand. 
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The  critical  Reynolds  number  for  laminar  instability  could  not  be  determined.  With  the 
assumption  of  rotational  symmetry  it  was  not  possible  to  confirm  the  experimental  observa- 
tion that  laminar  instability  occurs  at  a critical  Reynolds  number  of  approximately  130. 
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ABSTRACT 

A study  has  been  made  of  the  growth  through  collision  of  water  drops  in  the  atmosphere.  The  method 
of  superposition  of  flow  fields  obtained  from  the  numerical  solution  of  the  Navier-Stokes  equations  was 
used  for  the  calculations,  and  only  inertia,  gravity  and  drag  force  were  considered. 

The  calculated  linear  collision  efficiency  is  significantly  less  than  the  geometric  collision  efficiency 
Ik.- tween  a large  collector  drop  and  a small  collecting  drop  Itecause  of  the  effect  of  hydrodynamic  fortes 
The  linear  collision  efficiency  is  substantially  higher  than  the  geometric  collision  efficiency  between 
similarly  sized  drops  because  of  the  w-ake  effect. 

*1  o verify  the  validity  < if  the  calculations,  the  analytical  results  were  compared  with  available  ecperi 
mental  data.  Satisfa  tory  agreement  was  obtained  for  most  drop  sizes.  It  is  concluded  t hat  in  the  absence 
of  electricity  and  turbulence  the  dominant  factor  in  the  formation  of  precipitation  is  the  collisionat 
growth  of  similarly  sized  drops. 


1.  Introduction 

The  method  of  growth  of  water  drops  in  the  atmo- 
sphere can  he  classed  as  two  types:  condensational  anti 
collisional.  Condensational  growth,  in  which  water 
vapor  condenses  on  atmospheric  aerosols,  is  the  domi- 
nant mechanism  for  the  growth  of  water  drops  with 
radii  ft  20  pm.  The  W egener  Bergeron  Findeisen  pro- 
cess, in  which  the  ice  crystal  is  growing  at  the  expense  of 
supercooled  water  droplets  in  clouds  as  discussed  in 
standard  textbooks  such  as  Fleagleand  Businger  f 1 ‘>70) , 
is  also  considered  as  condensational  growth  of  ice  cry  - 
stals. Collisional  growth,  in  which  accretion  processes 
operate,  is  responsible  for  the  growth  of  water  drops 
with  radii  >20 pm.  This  study  is  concerned  with  the 
collisional  growth  of  water  drops,  because  it  is  directly 
related  to  the  understanding  of  the  precipitation  phe 
nomenon  in  clouds  as  well  as  to  the  scavenging  process 
of  pollutants. 

Analytical  studies  of  collisional  efficiency  have  been 
made  by  using  modi  lied  Stokes  flow  solutions,  modified 
• (seen  flow  solutions,  or  available  numerical  solutions  of 
the  Navier-Stokes  equations.  Owing  to  the  nonlinear 
nature  of  the  flow  held  surrounding  a relatively  large- 
drop,  which  is  responsible  for  the  collisional  growth, 
considerable  discrepancy  occurs  in  the  literature  for  the 
calculated  collision  efficiency.  Extensive  effort  has  been 
made  to  improve  the  calculation  procedures  for  collision 
efficiency  by  improving  the  empirical  relations  for  drag 
coefficients.  Based  on  an  improved  method  developed 
by  Lin  and  Lee  (1*27.1)  for  llow  field  calculations,  this 
study  shows  that  the  unsatisfactory  How  field  surround 


ing  a moving  water  drop  is  responsible  for  the  high 
order  of  magnitude  of  the  discrepancy  in  collision 
efficiency,  while  the  improvement  of  drag  formulas  can 
only  account  for  a very  small  percentage  of  the  dis 
crepancy  existing  in  the  literature. 

2.  Historical  background 

Collision  takes  place  when  one  moving  drop  (or 
collector  drop),  whose  movement  is  caused  In  gravita 
tional  or  other  forces,  collides  with  other  drops  (or 
collecting  drops)  which  happen  to  be  in  its  path.  To 
examine  the  probability  of  collision,  the  following 
equation  for  linear  collision  efficiency  is  used: 

J'r=yv  K ( 1) 

Here  K is  the  radius  of  the  collector  drop,  and  y,  is  the 
radius  of  the  collisional  cross  sectional  area  within 
which  any  collecting  drop  of  radius  r will  tollidc  with 
the  collector  drop.  Fig  1 shows  schematically  the 
nomenclature  used  in  this  study.  Ideally,  the  collisional 
cross-sectional  radius  should  be  equal  to  the  sum  of  the 
radii  of  the  collector  and  the  collecting  drops  if  the  dr  .ps 
are  noi  surrounded  by  air.  The  nondinicnsioiialized 
form  of  this  ideal  collision  cross  sectional  radius  is 
known  as  the  geometric  collision  efficiency.  However, 
as  a result  of  the  forces  generated  by  the  relative 
motion  between  the  water  drop  and  the  air,  laboratory 
observations  of  linear  collision  efficiency  often  differ 
substantially  yvith  the  geometric  collision  efficiency. 
These  phenomena  were  reported  by  T elford  ft  a /.  1 11). Tv, 
Schotland  and  Kaplin  (1*/56|,  Telford  and  T'h  undik 
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(1901).  Woods  (1965),  ami  Heard  and  I’ruppachar 
< 1‘WiX).  ID  understand  the  physical  processes  of  the 
collisionai  growth  of  water  drops,  many  investigators 
have  conducted  theoretical  studies  along  one  or  the 
other  of  the  following  two  approaches. 

a.  Evaluation  of  forcts  between  two  drops 

Hocking  (1959)  analyzed  the  Stokes  equations  and 
calculated  the  force  exerted  between  two  drops  to 
determine  collision  efliciency.  Davis  and  Sartor  (1%7) 
and  Hocking  and  Jonas  (1970)  extended  Hocking’s 
approach  to  include  rotation  of  the  water  drops.  How- 
ever, because  they  assumed  that  the  inertial  force  is 
negligibly  small  in  Stokes  equations,  their  approach  can 
only  apply  to  cases  of  practically  zero  Reynolds  num- 
ber ( Re~0).  For  freely  falling  water  drops  in  a standard 
atmosphere,  the  drop  sizes  and  terminal  falling  velocities 
vs  Reynolds  numbers  are  shown  in  Fig.  2.  It  can  be 
seen  that  a 30 /am  drop  has  a free  fall  velocity  of  10 
cm  s_l  that  corresponds  to  a Reynolds  number  of  0.4. 
The  effect  of  inertial  force  on  collision  efficiency  was 
considered  by  Klett  and  Davis  (1973)  using  a modified 
( tseen’s  equation  made  by  Carrier  (1953).  Their  results 
indicated  that  the  collision  efficiency  is  considerably 
larger  for  approximately  equal  size  drops.  However, 
due  to  the  inherital  limitation  of  the  Oseen’s  equation, 
this  method  is  only  valid  for  relatively  small  drops. 


h.  Superposition  of  flow  field 


Langmuir  (1948)  suggested  that  the  flow  field  around 
a single  drop  could  be  superposed  on  the  flow  held  of 
another  drop  to  determine  the  relative  motion  between 
the  two  drops.  Pearcey  and  Hill  (1956)  developed  this 
approach  by  using  ( loldstein’s  (1929)  solution  of  < Iseen’s 
(1910)  linearized  equations  for  flow  around  a spherical 
drop  at  a very  low  Reynolds  number  (Re<l).  Shafrir 
and  Neiburger  (1963),  Shafrir  (1965)  and  Xeiburger 
(1907)  improved  the  superposition  approach  by  using 
Jenson's  (1959)  method  to  obtain  solutions  of  the 
N’avier-Stokes  equations  for  an  intermediate  Reynolds 
number  range  (Re<20).  Shafrir  and  Tzvi  (1971)  ex- 
tended the  Reynolds  number  range  (Re<104)  by 
using  Rimon  and  Cheng’s  (1969)  method  with  a modi- 
fied boundary  condition.  Heard  and  drover  (1974) 


utilized  the  flow  field  i alculated  by  LcClair  et  at.  (1970), 
and  obtained  the  collision  efficiency  for  raindrops 
colliding  with  micron  size  particles.  Since  the  particle 
is  assumed  much  smaller  than  the  raindrops,  the 
„ existence  of  the  particle  does  not  alter  the  flow  held 

around  the  raindrop  (collector).  Consequently,  no 
superposition  is  necessary  if  the  radius  ratio  between 


20  /mi  collector  drop  and  a 19  /uni  collecting  drop  varies 
from  0 to  2.  To  resolve  some  of  these  discrepancies,  the 
available  literature  was  re-examined.  Hocking’s  ap 
proach,  which  appears  to  be  limited  to  very  small  drops, 
is  not  practical  when  applied  to  pollutant  scavenging 
and  weather  modification.  Langmuir’s  approach  does 
not  appear  to  have  any  size  restriction  for  most  atmo- 
spheric applications  provided  that  the  solution  of  the 
N’avier-Stokes  equations  can  be  accurately  determined. 
The  availability  of  high-speed  digital  computers  makes 
it  possible  both  scientifically  and  economically  to  meet 
this  requirement.  Le  Clair  et  at.  (1970)  and  Lin  and  Lee 
(1973)  obtained  numerical  solutions  for  flows  around 
spherical  drops  at  steady  and  transient  states,  respec- 
tively. The  accuracy  of  their  solutions  has  been  verified 
by  experimental  data  on  drag  coefficients,  separation 
angles,  and  wake  lengths.  This  study  uses  Langmuir’s 
approach  by  superpositioning  the  flow  fields  obtained 
from  Lin  and  Lee’s  (1973)  numerical  solution  of  the 
Navier-Stokes  equations. 

3.  Analysis 

Water  drops  in  the  atmosphere  may  be  influenced  by 
many  parameters,  such  as  gravity,  inertia,  turbulence, 
updraft  and  electrical  charge.  To  understand  the 
significance  of  each  parameter  on  droplet  growth, 
laboratory  experiments  have  to  be  conducted  in  a 
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controlled  environment.  Analytical  studies  are  being 
established  so  that  the  effect  on  droplet  growth  by  eat  it 
parameter  can  be  verified  with  experimental  data. 
With  a well-established  knowledge  concerning  each 
individual  parameter,  the  overall  effect  on  droplet 
growth  by  the  combination  of  the  many  physical 
parameters  actually  existing  in  the  atmosphere  tan 
then  be  investigated. 

If  one  considers  an  environment  where  only  gravity, 
inertia  and  drag  forces  are  significant,  the  equations  of 
motion  of  a freely  falling  drop  can  be  written  as 


hirfiR 

— (Vt-U.), 

l Ml 

6irgr 

- — (V.-Ui), 

. M. 


(2) 


(Ml 


in  which  the  subscripts  I.  and  .v  designate  the  larger 
collector  drop  and  the  smaller  collecting  drop,  tt-spet  - 
lively.  The  drop  velocity  is  V and  the  (low  field  \ < lot  it y 
U.  The  mass  of  the  water  tlrop  is  .V,  the  gravitational 
acceleration  is  g.  and  m is  the  dynamic  viscosity  The 
drag  coefficient  ( '/>  is  given  by  Oscen  (1910)  for  the 
Reynolds  number  region  of  0<Re$  1.6  as 

24 
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and  bv  Lin  and  Lee  (1973)  for  the  Reynolds  number 


A21 


C . 1. . I.  1 N ANI)  S . C . LK  K 


region  of  1.6<  Rc<  1000  as 
24 

Cn  — — (1+0.2207 


correction  formula 
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The  Reynolds  number  is  evaluated  by  using  the  relative 
velocity  between  the  drop  and  its  surrounding  flow 
held,  i.e., 

2p|VL-U.|R 

(Re)t=— , (6) 


2plV.-Ut|f 

t Re),= , 


in  wiii.  h p is  the  density  of  the  medium  in  which  the 
water  drop  is  traveling. 

To  calculate  the  position  of  each  drop,  Eqs.  (2) 
and  Id)  are  integrated  numerically  in  both  vertical  and 
horizontal  directions.  The  vertical  direction  is  along 
the  axis  connecting  the  center  of  the  collector  drop  and 
the  center  of  gravity  of  the  earth.  The  horizontal  direc- 
tion is  normal  to  the  vertical  direction  from  the  center 
of  the  collector  drop  toward  the  vertical  axis  of  the 
collecting  drop.  The  integration  is  performed  in  two 
steps  as  discussed  by  Conte  (1965).  The  first  step  is  an 
estimation  which  uses  the  Adams  Bashforth's  prediction 
formula 

V.H_v.++55(^Y-5++r 

24L  V ,tt  > \ ,11  / 

O"-*  1 

-9U  J- (8) 
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in  which  the  subscripts  n,  «+l,  ...,  denote  the  time 
steps  of  «,«+!,  ....  respectively.  The  velocity  at  the 
(>i+l)st  time  step  is  refined  by  the  Adams  . Moulton 


By  integrating  V with  respect  to  time  again,  the  position 
of  each  drop  can  be  determined.  These  processes  of 
integration  are  repeated  by  assuming  the  distance  of 
the  initial  separation  between  the  collector  and  the 
collecting  drops  of  given  sizes.  Collision  is  considered  to 
take  place  if  the  closest  distance  between  the  two  drops 
is  less  than  1%  of  the  smaller  collecting  drop  radius. 
An  initial  vertical  separation  distance  of  54  radii  of  the 
larger  collector  drop  is  used  because  the  flow  pattern 
calculations  indicate  that  the  tlow  fields  are  not  notice- 
ably disturbed  by  the  presence  of  each  other.  The  initial 
horizontal  separation  distance  is  the  variable  to  be 
determined.  The  largest  initial  horizontal  separation 
distance,  which  results  in  collision  between  the  collector 
and  the  collecting  drops,  is  known  as  the  radius  of  the 
coliisional  cross-sectional  area,  yc. 

4.  Results  and  discussion 

The  collision  trajectory  of  a collector  drop  and  a 
smaller  collecting  drop  is  shown  in  Fig.  4.  For  a collector 
drop  of  .40  gm  to  collide  with  a collecting  drop  of  9 gtn, 
the  coliisional  cross-sectional  radius  is  about  25  pm.  On 
the  other  hand,  the  40  /mi  collector  drop  can  only 
collide  with  a 4 gm  collecting  drop  within  a coliisional 
cross-sectional  radius  of  1.5  gm.  The  hydrodynamic 
force  in  the  vicinity  of  the  frontal  surface  of  the  collector 
drop  pushes  the  collecting  drop  away  from  its  falling 
axis.  The  collecting  drop,  as  a result  of  its  own  inertia, 
tends  to  maintain  its  own  course.  A larger  collecting 
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stantially  when  the  collector  drop  ^30  pm.  This  value 
differs  from  the  results  of  some  previous  investigators 
A comparison  of  theoretical  results  with  some  available 
experimental  data  is  shown  in  Fig.  7 for  the  case  of  a 
75  gm  collector  drop.  It  can  be  seen  that  the  agreement 
is  satisfactory  for  most  of  the  radius  ratio  regions  except 
for  those  drops  of  approximately  equal  size.  Re- 
examination of  the  theoretical  results  indicates  that  the 
vertical  distance  traveled  by  a pair  of  drops  of  appro.xi 
match  equal  size  is  very  long  before  the  collector  drop 
can  catch  up  with  the  collecting  drop.  Although  this 
particular  situation  exists  often  in  nature,  conventional 
laboratory  apparatus  is  limited  by  its  physical  dimen- 
sions and  prevents  this  phenomenon  from  being 
experimentally  observed. 

Comparison  is  also  made  with  the  analytical  results 
of  Klett  and  Davis  (1973)  for  a collector  drop  of  70  pm. 
It  is  noted  that  all  calculated  collisional  efficiencies  are 
large  for  high  radius  ratios.  However,  the  Klett  and 
Davis  results  appear  to  be  substantially  smaller  than 
available  experimental  data. 
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Fig.  5.  Collision  of  a pair  of  water  drops. 

drop  possesses  a greater  inertia  and  thereby  has  a 
larger  collisional  cross-sectional  radius. 

The  relative  position  of  two  drops  of  approximately 
equal  size  is  shown  in  Fig.  5.  For  a collector  drop  of 
194  pm  to  collide  with  a collecting  drop  of  183  (im,  the 
collisional  cross-sectional  radius  is  425  pm.  In  other 
words,  for  a pair  of  water  drops  with  a radius  ratio  of 
0.94.  the  linear  collision  efficiency  is  2.19,  which  is 
larger  than  the  geometric  collision  efficiency  of  1.94. 
The  reason  for  the  high  collision  efficiency  is  that  the 
disturbed  air  in  the  back  of  the  collecting  drop  drags  the 
collector  drop  into  its  wake  region.  The  disturbed  air 
in  the  front  of  the  collector  drop,  in  turn,  pushes  the 
collecting  drop  away  from  it.  [localise  the  horizontal 
distance  traveled  by  the  collector  drop  toward  the 
collecting  drop  is  larger  than  that  traveled  by  the 
collecting  drop  from  the  collector  drop,  the  pair  of 
drops  will  collide  even  though  their  initial  horizontal 
separation  is  larger  than  the  sum  of  their  radii.  This 
condition  is  sometimes  referred  to  as  the  wake  capture 
phenomenon,  w hich  is  one  of  the  im|xirtanl  processes  in 
prei  ipitation. 

The  variation  of  the  linear  collision  efficiency  with 
drop  radius  ratio  for  various  sizes  of  collector  drops  is 
shown  in  l ig.  f>.  It  should  be  noted  that  the  linear 
collision  efficiency  is  relatively  high  for  drops  of  approve 
mutely  equal  size.  The  numerical  value  increases  sub 
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Fig.  7.  Comparison  between  theory  and  experiment. 


5.  Conclusion 

The  linear  collision  efficiency  of  a pair  of  water  drops 
in  the  atmosphere  has  been  calculated  by  using  the 
method  of  superposition.  Based  on  the  (low  fields  ob 
tained  from  the  numerical  solution  of  the  Navier-Stokes 
equations,  collision  trajectories  betw  een  drops  of  various 
sizes  were  determined  by  considering  only  gravity, 
inertia  and  drag  forces.  For  cases  of  large  collector 
drops  and  small  collecting  drops,  the  linear  collision 
efficiency  is  significantly  less  than  the  geometric  collision 
efficiency,  because  the  hydrodynamic  force  upstream  of 
the  collector  drop  pushes  the  small  collecting  drops 
away  from  its  falling  axis.  For  cases  of  collector  and 
collecting  drops  of  similar  sizes,  the  linear  collision 
efficiency  is  substantially  higher  than  the  geometric 
collision  efficiency,  because  the  low  pressure  created  by 
the  viscous  forces  downstream  of  the  collecting  drop 
drags  the  collector  drop  into  its  wake.  The  theoretical 
result  has  been  compared  with  existing  data  from 
laboratory  experiments.  Satisfactory  agreement  be 
tween  the  theoretical  and  experimental  results  indicates 
that  the  theory  is  basically  sound.  In  an  atmospheric 
environment,  w here  only  inertia,  gravity  and  drag  forces 
are  significant,  the  dominant  factor  in  forming  precipita- 


tion appears  to  be  the  collisional  grow  th  of  drops  of 
water  of  approximately  equal  size.  Further  investiga- 
tions to  include  additional  atmospheric  parameters, 
such  as  electric  charges,  updraft  and  turbulence,  can 
be  developed  by  a similar  approach  for  applications  to 
more  complicated  atmospheric  problems. 
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APPENDIX 

Effect  of  Drag  Force  on  Collision  Efficiency 

A water  drop  in  the  atmosphere  is  surrounded  by  air, 
which  produces  a drag  force  when  the  drop  travels  w ith 
a certain  velocity.  Calculation  of  collision  efficiency 
requires  a knowledge  of  the  drag  force  of  various  drop 
sizes  at  different  velocities.  Consequently,  a dimension 
less  drag  coefficient  is  given  as  a function  of  Reynolds 
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numbers.  At  a Reynolds  number  region  between  0 and  1, 
it  is  not  absolutely  certain  that  the  drag  coefficient 
actually  follows  the  Stokes  drag  formula  or  Oseen’s 
drag  formula.  To  examine  the  actual  derivation  of  the 
linear  collision  efficiency  caused  by  the  uncertainty  of 
the  chosen  drag  relation,  Fig.  8 shows  the  comparison 
of  the  Stokes  and  Oseen’s  drag  on  linear  collision 
efficiency  calculations  for  a 30 /*m  collector  drop.  It  is 
obvious  that  the  drag  formula  does  not  affect  the  trend 
of  the  linear  collision  efficiency.  A maximum  deviation 
of  10%  could  be  expected  in  the  vicinity  of  drop  radius 
ratio  of  0.75.  For  larger  drops,  evaluations  of  collision 
efficiencies  were  made  by  using  Heard  and  Fruppacher’s 
formulas  as  well  as  Lin  and  Lee’s  formula  based  on  the 
same  flow  held  data  for  a drop  radius  of  100 /un.  1 ig.  8 
shows  that  no  appreciable  difference  can  be  found  in 
linear  collision  efficiency  for  radius  ratios  <0.7.  In  the 
region  where  the  sizes  of  collecting  drops  are  approach- 
ing the  size  of  the  collector  drop,  no  drag  formula  can 
cause  the  decrease  of  collision  efficiency  as  reported  in 
some  early  investigations.  Step  changes  in  the  drag 
formula  are  the  only  cause  of  the  irregularities  of  the 
collision  efficiency  curve,  and  these  are  not  readily 
explainable  in  natural  phenomena. 
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In  response  to  Dr.  Klett’s  comment  on  the  article  of 
Lin  and  Lee  (1975),  it  is  necessary  to  keep  the  following 
subjects  in  focus. 

1.  The  purpose  of  a scientific  article 

A scientific  article  is  being  published  for  the  purpose 
of  clarifying  a situation  which  is  being  confused  by 
confiicting  results,  bigs.  1 and  2 show  the  situation  on 
linear  collision  efficiencies  today  for  a collector  drop  of 
approximately  70  ^m  and  30 /im,  respectively.  It  is 
evident  that  the  method  of  superposition  used  by  Lin 
and  Lee  (1975)  and  the  method  of  force  evaluation  used 
by  Klett  and  Davis  (1973)  provide  better  agreement 
in  comparison  with  the  results  in  earlier  studies.  To 
argue  for  the  better  accuracy  of  one  method  over  the 
other  does  not  serve  any  constructive  purpose,  because 
both  are  approximate  in  nature.  [This  has  been  pointed 


out  by  Lin  and  Lee  (1975)  as  well  as  by  the  original 
authors  of  these  methods,  Langmuir  (1948)  and 
Hocking  (1959).]  Improvement  by  Lin  and  Lee  (1973) 
in  computation  techniques  can  only  minimize  the 
inaccuracy  which  is  caused  by  the  inability  of  obtaining 
mathematical  solutions  of  the  Navier-Stokes  equations 
applying  to  a flow  f.jld  that  consists  of  two  or  more 
drops  of  various  sizes.  Similarly,  modification  of  Oseen’s 
linearization  technique  by  compromising  with  Stokes’ 
approach,  as  discussed  by  Carrier  (1953)  and  Klett  and 
Davis  (1973),  is  not  a solution  of  the  nonlinear  Navier- 
Stokes  equations  which  govern  flow  fields  of  continuum 
media  at  non-zero  Reynolds  numbers.  A 70  jim  radius 
water  drop  falling  in  air  at  steady  state  gives  a Reynolds 
number  of  approximately  4,  and  a 30  /im  drop  of 
approximately  0.4.  When  the  flow  field  is  being  dis- 
turbed by  another  drop  which  occurs  before  all  colli- 
sions, the  assumption  of  steady-state  flow  phenomenon 
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Fig  1 Linear  collision  efficiency  for  collector  drop 
radius  of  approximately  70  Mm. 


FlC.  2.  Linear  collision  efficiency  for  collector  drop 
radius  of  approximately  30  Mm. 
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is  in  question,  not  to  mention  the  credibility  of  super- 
position or  linearization.  It  would  be  a worse  alternative 
if  this  discussion  were  to  leave  the  impression  that  the 
problems  caused  by  superposition  could  be  eliminated 
by  some  kind  of  compromise  between  two  linearization 
techniques. 

2.  A curve  on  linear  collision  efficiency 

In  Fig.  7 of  Lin  and  Lee  (1975),  a curve  on  collision 
efficiency  E by  Klett  and  Davis  (1975)  was  used  as 
linear  collision  efficiency  Yc  for  comparison  between 
theory  and  experiment.  The  comparison  showed  that 
the  numerical  values  in  the  analytical  result  of  Klett 
and  Davis  were  less  than  those  of  Lin  and  Lee,  which 
agreed  well  with  available  experimental  data.  This  was 
a mistake  by  Lin  and  Lee.  The  result  of  Klett  and  Davis 
was  presented  as  collision  efficiency  not  as  linear 
collision  efficiency.  A sincere  apology  is  extended  to 
Klett  and  Davis  for  the  hasty  addition  of  that  particular 
curve.  Our  explanation  here  is  not  to  find  an  excuse 
about  what  was  done  but  to  present  the  circumstances 
of  turning  a good  intention  to  a poor  relation.  The 
manuscript  of  Lin  and  Lee  was  prepared  before  the 
article  of  Klett  and  Davis  appeared  in  the  Journal.  The 
question  in  doubt,  at  that  time,  was  “whether  the  wake 
capture  phenomenon  increases  or  decreases  collision 
efficiency.”  Most  earlier  publications  (Davis  and 
Sartor,  1967;  Shafrir  and  Neiburger,  1963;  Hocking 
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Fjg.  4.  Comparison  of  collision  efficiencies  obtained  b\ 
superposition  and  force  methods. 


and  Jonas,  1970)  indicated  that  the  collision  efficiencies 
tended  to  go  to  zero  as  the  droplet  radius  ratios  ap- 
proached unity.  Lin  and  Lee  found  that  the  trend  was 
exactly  the  opposite.  The  trend  of  Klett  and  Davis 
agreed  with  that  of  Lin  and  Lee  and  was  quickly  added 
to  Fig.  7 during  che  revision  of  their  article  for  sup- 
porting the  controversy.  It  was  not  conceivable  that 
those  two  curves,  which  were  calculated  from  two 
completely  different  methods,  could  be  in  agreement 
within  a few1  percent  at  a time  when  earlier  studies  had 
discrepancies  ranged  almost  between  zero  and  infinity. 
In  order  to  clarify  the  situation  on  collision  efficiency 
analyses,  the  positive  aspect  of  these  two  artioes  needs 
to  be  emphasized  and  is  shown  in  Figs.  3 and  4 for  the 
linear  collision  efficiency  and  the  collision  efficiency 
respectively.  The  methods  of  superposition  and  force 
evaluation  are  approximations  which  produce  com 
parable  results  for  collision  efficiency  analyses. 

3.  Oseen’s  drag  equation 

The  Oseen’s  drag  equation  was  given  as  Eq.  (4)  by 
Lin  and  Lee  (1975).  A table  inserted  in  Fig  8 of  the 
article  misprinted  the  same  equation  by  omitting  a 
Reynolds  number.  This  was  an  oversight  by  1 in  and 
Lee.  The  thoroughness  of  Dr.  Klett  in  reading  this 
article  is  to  be  commended.  We  appreciate  his  prompt 
ness  in  bringing  it  to  our  attention. 
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The  experimental  data  of  Cataneo  rt  al.  (1071)  and 
the  analytical  results  of  Lin  am!  Lee  (1175)  appear  to 
be  in  excellent  agreement  on  collision  efficiency,  which 

1 Present  affiliation  National  Aviation  I'acilitics  K\|it-rinicntal 
Center.  Atlantic  City,  N J 08405 


is  the  nondimensionali/.ed  distance  of  the  maximum 
horizontal  separation  between  two  free-falling  drops. 
However,  an  apparent  discrepance  exists  with  regard 
to  the  amount  of  their  vertical  separation  that  max 
affect  the  collision  efficient-;. . 1 pon  reexamining  the  two 
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papers,  the  following  distinctions  are  noted  : 

1)  The  study  of  Lin  and  Lee  (1975)  was  undertaken 
to  determine  collision  efficiency  by  using  a vertical 
separation  distance,  which  is  critical  for  the  accuracy 
of  tire  calculation.  A maximum  horizontal  separation 
distance,  which  w ill  result  in  the  collision  of  two  falling 
drops  of  approximately  ecjua!  size,  occurs  within  the 
range  of  54  radii  of  the  vertical  separation  distance.  In 
other  words,  the  etfect  on  collision  efficiency  calculations 
becomes  insignificant  when  the  vertical  separation 
distance  is  largt  - than  54  radii.  It  is  necessary  to  point 
out  that  this  hypothesis  was  only  indirectly  ex, duated 
through  the  effect  on  (he  drag  coefficient  by  assuming 
various  sizes  of  the  outer  boundary,  which  approximates 
the  location  of  the  undisturbed  tlow  field  in  numerical 
calculations.  Detailed  discussions  of  the  effect  of  the 
outer  boundary  on  drag  coefficient  calculations  have 
been  published  b\  Hamielec  cl  a/.  (1967),  Let  lair  et  ol. 
(1970)  and  Lin  and  Lee  (1975).  Fig.  1 illustrates  the 
test  range  of  Cataneo  et  ill.  (1971).  Reynolds  numbers 
Re  of  2(1  and  40  are  approximately  equal  to  free-falling 
drops  w ill)  radii  of  140  and  ISO  gm,  respet  lively  . It  r an 
he  seen  that  the  assumption  of  54  radii  appears  to  he  a 
conservative  estimate  for  drag  coefficient  evaluations. 
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2)  In  the  experiment  of  Cataneo  el  <il.  (1071)  for 
drops  of  equal  si/e,  the  wake  el  feet  of  the  collecting  drop 
was  observed  at  a vertical  separation  distance  > KM) 
diameters  because  the  collector  drop  was  noticeably 
;u t derate 1 1.  In  order  t > compare  t heir  observation  with 
the  theorel i<  al  value  of  54  radii  that  is  used  in  collision 
efficients  calculations,  it  is  necessary  to  ask  whether 
the  maximum  horizontal  separation  distance  is  smaller 
or  larger  than  the  critical  limit  determined  from  the 
collision  elite  ic-nt  \ . In  other  words,  although  the  wake 
length  m.is  lu*  longer  than  KM)  diameters,  the  elfee  t on 
the  calculated  collision  elite  iem  \ Incomes  insignificant, 
as  long  as  the  radius  of  the  wake  is  less  than  the  maxi- 
mum horizontal  separation  distance  that  determines 
the  collision  efficiency . 

The  excellent  agreement  on  collision  elite  iem  s be- 
tween the  experiment  of  ( alamo  el  al.  (ULl)  and  the 
thcorv  of  Lin  and  Lee  (1^75)  indicates  that  the  hy- 
pothesis is  reasonable.  Direct  verification  of  this 
hypothesis  is  desirable;  however,  a numerical  investiga- 
tion of  the  vertical  separation  is  impractical  because 
too  much  computer  time  is  required  for  calculating 
slreamfunc  lions  and  vorticitics  of  an  extraordinarily 
large  flow  field  and  for  evaluating  by  the  method  of 
superposition  the  change  of  vertical  separation  from 
two  seemingly  undisturbed  flow  lields.  h xperimental 
investigations  may  be-  more  realistic  if  existing  facilities 
are  used  lo  measure  photographically  the  horizontal 
separation  simultaneously  with  the  vertical  separation. 
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HEAT  MASS  AND  MOMENTUM  TRANSFER  IN  TURBULENT  BOUNDARY  LAYER  FLOWS 

S.C.  Lee,  D.W.  Pepper,  W.M.  Byrne,  and  R.C.  Tai 
University  oi  Missouri,  Rolla,  U.S.A. 


Abstract 

An  analytical  method  using  turbulence  energy  to  study  the  transfer  of  heat  mass  and  momentum  in  a turbu- 
lent boundary  layer  was  developed.  This  approach  considers  the  history  of  the  turbulent  motion  by  treat- 
ing the  eddy  viscosity  as  a dependent  variable  to  be  determined  as  a function  of  locations  in  a turbulent 
flow  field.  A systematic  comparison  with  experimental  data  was  made.  The  same  empirical  constants  es- 
tablished in  an  incompressible  self-preserved  turbulent  flow  along  a flat  plate  were  also  applicable  tc 
a:ceierat  ing  and  decelerating  flows  with  suction  and  blowing  as  well  as  tor  boundary  layer  flows  with 
heat  and  mass  transfer. 


NOMENCLATURE 


C , c 
D 
F 
g 

H , h 
K 
P 
P 

P+ 

Pr 

Re 

Ri 

St 

T 

U,  u 

V* 

V,  V 


> 

a' s 
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Species  concentration*  rag^/cc  ^ 
Turbulence  dissipation,  m“/sec 
Dimensionless  blowing  parameter 


2 


Gravitational  acceleration,  m/sec 
Enthalpy,  nr /sec z 2 

Kinetic  energy  of  turbulence,  m /sec 
Turbulence  production,  m^/sec^ 

Static  pressure,  newton/tn* 

Dimensionless  pressure  gradient 

Prandtl  number 

Reynolds  number 

Richardson  number 

Stanton  number 

Temperature,  °k 

Velocity  component  in  /.-direction,  m/sec 

Friction  velocity,  m/sec 

Velocity  component  in  y-direction,  m/sec 

Horizontal  distance,  m 

Vertical  distance,  m 

Dimensionless  vertical  distance 

Coefficients  of  the  generalized  equation 

Coefficient  in  the  modified  law  of  the  wall 

Boundary  layer  thickness 

Exchange'  coefficient,  newton-sec  /m^ 

Von  Karman’s  universal  constant 

2 

Dynamic  viscosity,  newton- s ec /m 
Kinematic  viscosity,  m /sec 
Density,  newton-sec^ /m^* 

Dimensionless  constant., 

Shear  stress,  newton/m"’ 

Generalized  flow  parameter  2 
Stream  function,  newton-sec /ra 
Normalized  stream  function 


Subscr Ipt 

C : Concentration 

H : Energy 

K Kinetic  energy  of  turbulence 

M : Momentum 

U : Wall 

0 : Initial  x location 

00  : Free  stream 


INTRODUCTION 

One  of  the  most  commonly  observed  phenomena  occur- 
ring in  many  practical  problems  is  the  transfer  oC 
heat,  mass,  and  momentum.  Most  frequently  this 
phenomenon  is  encountered  in  turbulent  flows  w .ich 
are  governed  by  the  Reynolds  Equations.  Analyti- 


cal solutions  of  the  Reynolds  equations  art  only 
tb  retically  possible  if  the  turbulent  stresses 
known.  In  order  to  close  the  Revnolds  equa- 
.ons,  Boussinesq  / 1 / related  the  turbulent  shear 
stress  with  the  average  velocity  gradient  bv  intt 
ducing  an  "eddy  viscosity,"  also  known  as  a ’momen- 
tum exchange  coefficient."  It  is  the  precise  l.*rn 
of  this  "eddy  viscosity"  that  becomes  the  sub'*  t 
of  many  phenomena  1 ogica 1 theories.  Prandtl  121 
used  momentum  as  a transferable  quantity  resulting 
in  a mixing  length  which  was  considered  as  the  av- 
erage distance  between  turbulent  eddies  in  analogy 
to  the  mean  free  path  between  molecules  in  kinetic 
theory.  Taylor  / 3/  used  the  vorticity  as  a trans- 
ferable quantity  and  also  derived  a similar  length 
parameter.  Nevertheless,  neither  Prandtl  nor  Tav- 
lor  could  give  a workable  relation  for  this  length 
parameter  because  of  the  lark  of  a universal  con- 
stant. Von  Karman  /<*/  conducted  a series  oi  ex- 
periments through  which  a length  parameter  was 
found  to  be  linearly  proportional  to  the  ratio  of 
the  first  and  second  derivatives  of  the  average 
velocity.  He  named  his  value  of  proportionality 
the  "univ.-rsal  constant"  which  was  applicable  to 
some  turbu’ent  flow  problems.  Earlier  eifcrts  in 
establishing  a workable  phenomenological  theory  oi 
turbulence  were  discussed  in  detail  by  Hinze  / 5/ - 
It  is  evident  that  the  mixing  length  approach  did 
not  provide  the  necessary  mechanism  to  take  the 
history  of  turbulence  into  account.  The  use  t 
the  turbulence  kinetic  energy  equation  to  complete 
the  Reynolds  equations  was  first  proposed  by  Kol- 
mogorov /6/  and  discussed  In  standard  textbooks, 
such  as  Townsend  111  % Hinze  / 5 / , Monin  and  Yaglom 
/ 8/  and  many  others.  Kolmogorov  / 6 / and  Prandtl 
/ 9/  independently  suggested  that  the  turbulent 
shear  stress  related  to  the  turbulence  kinetic  en- 
ergy's one-half  power.  This  model  has  been  the 
basis  of  calculation  procedures  by  Glushko  110/ 
and  Spalding  / 1 1 / . A simpler  model  was  suggested 
by  Nevzglajdov  / 1 2 / in  which  the  turbulent  shear 
stress  related  linearly  to  the  turbulence  kinetic 
energy.  The  latter  model  was  used  by  Bradshaw  et 
al  / 1 3 / and  Lee  and  Harsha  /\U/  for  boundarv  layer 
flows  and  free  turbulent  flows,  respectively. 

Figure  1 shows  the  correlation  of  the  Nevzglajdov 
model  with  the  experimental  data  of  Klebanoft  ,15/ 
and  Schon  and  Mery  / 1 6 / . It  is  evident  that  the 
linear  relation  between  turbulent  shear  stress  and 
turbulence  kinetic  energy  describes  realistically 
the  physical  phenomenon.  Lee  et  al  / 1 7 / applied 
the  Nevzglajdov  model  successfully  in  studving 
heat  mass  and  momentum  transfer  in  tree  nixing 
flows.  This  paper  is  to  extend  the  same  approach 
to  boundary  layer  flows. 
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Fig.  1:  Data  Correlation  of  Nevzglajdov  Model 
ANALYSIS 

The  mathematical  formulation  of  heat,  mass,  and  mo- 
mentum transfer  in  turbulent  boundary  layer  flows 
was  given  by  Patankar  and  Spalding  / 1 8 / . The  use 
of  Nevzglajdov  model  to  relate  the  turbulent  shear 
stress  and  the  turbulence  kinetic  energy  was  dis- 
cussed by  Lee  and  Harsha  / 1 4 / . The  validity  of 
using  constant  turbulent  Prandtl  and  Schmidt  num- 
bers was  examined  by  Lee  et  al  / 1 7 / . A brief  sum- 
mary of  the  analytical  approach  of  the  turbulence 
kinetic  energy  method  for  various  boundary  layer 
flow  problems  can  be  outlined  as  follows: 

Governing  Equations 

For  two-dimensional,  steady  state,  boundary  layer 
flows,  the  continuity,  momentum,  turbulence  kine- 
tic energy,  mean-flow  energy  and  species  concentra- 
t ion  equations  may  be  written  is: 
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where  U and  V are  the  time-average  velocity  com- 


ponents in  the  x and  y directions,  respectively;  p 
is  the  density;  p is  the  static  pressure,  K is  the 
kinetic  energy  of  turbulence;  H is  the  stagnation 
enthalpy;  and  C is  the  species  concentration.  Us- 
ing the  Nevzglajdov  model  as  modified  by  Lee  and 
Harsha  / 1 4 / , the  turbulent  shear  stress,  T,  can  be 
related  to  the  turbulence  kinetic  energy,  K,  as 
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The  exchange  coefficients  may  then  be  defined  as: 
Momentum : 

eM--<pvU>/g-0.3pK/||H|  (7) 

Turbulence  Kinetic  Energy: 

3K 

CK  - -<pvK  > I-  . rw/oK  (8) 

Mean-Flow  Energy: 

EH  * '<PVh  ” ^ ' fM/0H  (9> 

Species  Concentration: 

°c  * '<PVC  > - cM/oc  (10) 

where  the  lower  letter  case  designates  the  fluctua- 
ting components.  Og  is  the  equivalent  Prandtl  num- 
ber which  is  approximately  equal  to  0.7  determined 
by  Lee  et  al  /I 7/.  is  the  turbulent  Prandtl 

number  with  a numerical  value  of  0.75  obtained  ex- 
perimentally by  Arya  / 1 9 / . o^.  is  the  turbulent 

Schmidt  number  with  a numerical  value  of  0.75  given 
by  Fleagle  and  Bussinger  / 2 0 / for  air  and  water  va- 
por mixture  in  the  atmosphere.  The  turbulence  pro- 
duction term,  P^,,  consists  of  the  shear  generated 
turbulence  due  to  the  mean  motion  and  the  buoyancy 
generated  turbulence  due  to  flow  stratification. 

As  used  by  Plate  / 2 1 / , the  production  term  may  be 
written  as 
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with  g designating  the  gravitational  acceleration 
and  Ri  designating  the  Richardson  number.  The  dis- 
sipation term,  D^.,  was  first  used  by  Glushko  / 1 0 / 
and  later  modified  by  Byrne  1221  as 
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3/2 

Dr  - 1.8pKJ/  /6 

3/2 

Dr  - (<5/4y)  1 .8pKJ/  /6 


for  y > 5/4 
for  y < 5/4 


(12) 


(13) 


where  6 is  the  boundary  laver  thickness.  Bvrne's 
expression  is  in  close  agreement  with  that  of  Brad- 
shaw et  al  / 1 3 / for  the  region  near  the  wall  boun- 
dary . 

Numerical  Method 

Except  the  continuity,  all  governing  differential 
equations  are  parabolic.  Transforming  to  stream- 
line coordinates,  (x,ii>),  the  system  of  equation  can 
be  solved  with  the  numerical  method  developed  bv 
Spalding  and  Patankar  / 2 3 / modified  by  Lee  and  Har- 
sha / 14  / . The  generalized  parabolic  equation  takes 
the  form: 
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Table  1:  Coefficient  of  of  Eq.  (14) 
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where  u>  is  the  normalized  str eamf unc t ion  with  to  e 
qual  to  zero  or  unity  at  the  wall  boundary,  «^w,  o 
the  free  stream  boundary,  , respectively.  ihe 
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and  can  be  expressed  as 
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where  0A  » 1,  o , 0^  or  oc  for  <t>  ■ U,  K,  H,  or  C, 
respectively,  the  coefficient  for  various  val- 
ues of  4>  is  given  in  Table  1. 
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Fig.  3:  Heat  Transfer  zn  Turb.  Boundary  Layer 
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Boundary  Cund 1 t Ions 

Standard  boundary  conditions  are  being  used  at  the 
free  stream  boundary  as  well  as  at  the  wall  bound- 
ary. However,  in  the  laminar  sublayer  adjacent  to 


V ^/U*  P*  eq . (1C)  exp.  ref. 
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the  wall  boundary,  the  turbulent  shear  stress  de- 
creases rapidly  to  zero  while  the  laminar  shear 
stress  increases  to  the  value  of  the  wall  shear, 
Ty.  Defining  an  "effective  shear  stress"  as  the 
summation  of  the  local  laminar  and  turbulent  shear 
stresses,  an  "effective  turbulence  kinetic  energy" 
can  be  introduced  to  satisfy  the  linear  relation 
of  the  Nevzglajdov  model.  Using  the  "slip  value" 
approach  of  Spalding  and  Patankar  / 23,',  the  fric- 
tion velocity,  U* , can  be  related  to  the  local  ve- 
locity, U,  by  a modified  "Law  of  the  Wall": 


jjj  - 12-5  lny+  + B][l  + (2.5  In/  + 6)1  (16) 
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correlation  of  Eq.  (16)  with  the  experimental  data 
of  Rotta  / 2 4 / and  Julien  et  al  / 25/  is  shown  in 
Figure  2 for  boundary  layer  blowing  and  suction 
under  the  influence  of  pressure  gradient,  \ 
Consideration  of  thermal  stratifications  was  bast 
on  the  suggestion  of  Lumley  and  Panofskv  12b/, 
where  Ri  is  the  Richardson  number.  In  heat  tr.ms 
fer  studies,  either  the  Nusselt  number  >r  the 
Stanton  number  is  being  used.  Modifying  the  *x 
pression  of  Reichardt  / 2 7 / , the  Stanton  T.-.n.hc.  , 

St,  becomes: 
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with  F ■=  (pV)u/(pU)tIJ 

For  Che  case  of  Py  equal  to  0.75  and  Pr  equal  to 
0.72,  a is  approximately  equal  to  0.96.  Correla- 
tion of  Eq.  (17)  with  the  experimental  data  of 
Moffat  and  Kays  / 28/  is  shown  in  Figure  3. 
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RESULTS  AND  DISCUSSION 


Application  of  Che  turbulence  kinetic  energy  meth- 
od was  made  syxtenat ica 1 ly  to  a number  of  engi- 
neering problems.  lai  129/  analyzed  turbulent 
boundary  layer  flows  along  a flat  plate  and  com- 
pared with  the*  experimental  data  of  Klebanoff  / 15/. 
Byrne  1221  analyzed  turbulent  boundary  layer  flows 
with  blowing  and  suction  under  the  influence  of 
favorable  and  adverse  pressure  gradients  and  com- 
pared with  the  experimental  data  of  Julien  / 30 / 
and  Bradshaw  / 31/.  Figure  U shows  one  of  Byrne's 
comparisons  for  boundary  layer  blowing  with  fa- 
vorable pressure  gradient.  Good  agreement  was  ob- 
tained between  Byrne's  analytical  results  and 
Jullen's  experimental  data  in  velocity  profiles. 
Only  sheer  stress  at  the  wall  was  measured.  As 
expected,  the  maximum  sheer  stress  for  the  blowing 
case  occurred  at  a small  distance  avay  from  the 
wall.  Pepper  / J 2 / analyzed  thermally  stratified 
boundary  layer  in  the  atmosphere  and  compared  his 
results  with  some  available  wind  tunnel  data.  Fig- 
ure 5 shows  the  comparison  of  Pepper's  analytical 
results  with  Arya's  719/  experimental  data  on  a 


: Velocity, 

Temperature , 

and  Turb. 

Enerqy 

x/6 

o 

theory 

exp. 

ref. 

13.1 

o 

19 

16.6 

— 

A 

19 

stable  atmosphere  boundary  layer  (cold  surfa*  . 
Good  agreement  was  obtained  for  velocity,  tempera- 
ture and  turbulence  energv  profiles.  Figure  b 
shows  the  comparison  of  Pepper's  results  with  Nal-V 
hotra’s  / 3 3 / experimental  data  on  a simulated  pol- 
lutant (ammonia)  diffusion  process  in  an  unstable 
atmosphere.  Good  agreement  was  obtained  for  veloc- 
ity, temperature  and  concent  rat  ion  d ist r i but  ions . 

CONCLUSION 

The  turbulence  kinetic  energv  equation  was  used  to 
complete  the  Reynolds  equations  for  studying  the 
heat  mass  and  momentum  transfer  in  turbulent 
boundary  layer  flows.  By  calculating  the  eddy 
viscosity  everywhere  in  the  flow  field,  a set  of 
uniquely  determined  empirical  constants  was  estab- 
lished for  boundary  laver  1 lows  along  a tlat  plate 
with  zero  pressure  gradient,  for  accelerating  and 
decelerating  flows  with  suction  and  blowing,  as 
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well  as  for  heating  and  cooling  of  the  wall  bound- 
aries. This  approach  gives  realistic  solutions 
not  only  to  the  engineering  design  parameters  of 
wall  shear  and  surface  heat  transfer  but  also  to  the 
the  meteorological  parameters  which  may  be  used  for 
studying  pollutant  diffusion  in  the  atmospheric 
boundary  layer. 
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Transport  Phenomena  in  Thermally 
Stratified  Boundary  Layers 

Studies  of  heat,  mass,  and  momentum  transfer  are  made  to  analyze  the  formation  of  ma- 
rine fogs  in  thermally  stratified  boundary  layers  in  the  atmosphere  I he  governing  par 
tial  differential  equations  of  continuity,  momentum,  temperature,  and  concentration  are 
used  to  describe  the  transport  phenomena  An  additional  equation  of  turbulence  energy 
is  introduced  to  account  for  the  development  of  the  turbulent  motions  Simultaneous  so- 
lution of  this  system  of  equations  allows  the  turbulent  exchange  coefficients  to  be  treat 
ed  in  the  same  way  as  all  other  dependent  parameters  Verification  of  the  theoretical 
approach  is  made  by  comparing  the  numerical  predictions  with  wind  tunnel  simulations 
of  neutral,  stable,  and  unstable  atmospheres  Application  of  the  theory  is  extended  to 
the  investigation  of  the  formation  of  advection  fog  over  cold  ocean  surfaces  In  addition 
to  the  established  criteria  obtained  from  wind  tunnel  data,  the  fog  model  takes  into  con 
sideration  the  radiation  and  sedimentation  of  fog  droplets  as  veil  as  condensat ion  and 
evaporation  of  liquid  water 


1 Introduction 

Considerable  effort  has  been  made  in  recent  years  to  understand 
transport  phenomena  in  thermally  stratified  boundary  layers  Ad 
equate  knowledge  in  this  area  is  valuable  for  predicting  the  diflu 
sion  process  of  air  pollutants  in  the  lower  atmosphere  as  well  as  for 
forecasting  air  water  circulation  for  weather  conditions.  Manv 
physical  parameters  are  involved  in  atmospheric  transport  pro 
cesses,  such  as  wind,  temperature,  and  concentration  of  the  diffus 
ing  medium  as  well  as  the  geographical  terrain  Studies  of  at  mo 
spheric  motions,  however,  are  hindered  by  the  turbulence  generat 
cd  from  the  interactions  of  all  these  related  parameters.  Moreover, 
as  a result  of  the  random  motion  of  turbulence  eddies,  field  mea 
surementa  are  often  insufficient  for  formulating  a reliable  mat  he 
matical  model 

Simulations  of  atmospheric  boundary  layers  are  usually  accom 
plished  by  wind  tunnel  modeling  and  numerical  computations.  A 
detailed  survey  of  the  literature  concerning  field  measurements, 
laboratory  simulations,  and  numerical  predictions  is  given  by  Pep 
per  |l|  Although  many  questions  have  been  raised  as  to  whether 
atmospheric  turbulence  can  be  realistically  simulated  in  a wind 
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tunnel,  no  definite  answer  is  readily  available  As  a consequence,  it 
appears  necessary  to  develop  a realistic  model  to  numerically  ana 
lyze  transport  phenomena  which  not  only  can  he  verified  bv  wind 
tunnel  experiments  but  may  also  be  used  for  predicting  at  mo 
spheric  motions. 

Mathematical  analyses  of  turbulent  transport  processes  are  dp- 
cussed  by  Hinze  (2),  Monin  and  Yaglom  |M| . and  Lumley  and  Pa 
nofsky  f 4 1 Based  upon  their  studies,  transport  phenomena  are 
normally  treated  by  phenomenological  and  statistical  theories 
While  physically  more  realistic  than  phenomenological  models,  the 
use  of  statistical  description  requires  a large  volume  of  carefully 
sampled  experimental  data  which  are  not  currently  available.  Phe 
nomenological  theory  is  based  principally  upon  mixing  length 
theories  and  turbulence  kinetic  energy  approaches.  The  mixing 
length  models,  based  upon  the  concept  developed  by  Prandtl  |ft|, 
incorporate  rather  simple  yet  successful  empirical  relations  to  av 
count  for  the  Reynolds  shear  stress  term  in  the  equations  ot  mo 
tion.  Phe  use  of  turbulence  kinetic  energy  to  complete  the  Reyn 
olds  equations,  first  proposed  by  Kolmogorov  |6|.  is  discussed  in 
textbooks  by  Townsend  |7],  Hinze  |2|.  Monin  and  Yaglom  \W\. 
Tennekes  ano  Lumley  [8],  and  many  others  Following  the  sugges 
Hons  of  Kolmogorov  [fi|  and  Prandtl  |f»|.  the  turbulent  shear  stress 
is  related  to  the  turbulence  kinetic  energy  's  one-hall  power  I his 
model  is  the  basis  of  calculational  procedures  used  by  (ilushko  ' 
Patankar  and  Spalding  |10|.  and  Goaman.  et  a)  |JJ|  A simple -r  re 
lationship  is  suggested  bv  Nevzglajdov  [ 12 1 in  which  the  turbulent 
shear  stress  is  linearly  related  t«  the  turbulence  kinetu  energy 
This  concept  has  been  successfully  used  bv  Bradshaw,  el  al 
and  Lee  and  Harsha  1 1 4)  in  analyzing  boundary  layer  and  free  tur 
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hulent  flows,  respectively 

In  lieu  of  the  many  different  approaches  used  in  t losing  the  dll 
ferential  equations  governing  turbulent  flows,  it  is  necessary  to  in 
vestigate  the  validity  of  all  available  models  before  conducting  an 
analysis  Bv  using  the  experimental  data  of  Klebanoff  ( 1 f> | and  the 
neutral  wind  tunnel  data  of  Schon  and  Mery  J lb)  as  a basis  of  com 
parison.  it  was  possible  to  analyze  the  three  methods  of  closure 
that  are  shown  in  Fig  1 The  disarray  of  data  points  in  the  Prandtl 
mixing  length  model  indicates  that  the  mixing  length  parameter 
requires  modification  for  each  particular  problem.  The  Kolmogo 
rov  model,  while  more  adequate  than  the  Prandtl  mixing  length 
model.  Iiegins  to  show  some  scattering  in  the  near  wall  regions  of 
the  boundary  layer  The  Nevzglajdov  model,  on  the  other  hand,  is 
simple  and  gives  relatively  less  scattering  lor  the  same  data.  Sirni 
lar  conclusions  are  also  drawn  by  Harsha  and  Lee  |17j  for  free  tur- 
bulent flows. 

The  approach  adopted  in  this  study  attempts  to  make  optimum 
use  of  the  current  advancement  in  numerical  techniques  coupled 
with  presently  available  experimental  data  for  the  purpose  of  real 
istically  predicting  atmospheric  motions  with  a minimum  of  em 
piricism.  Comparison  of  the  numerical  results  with  wind  tunnel 
simulations  under  various  conditions  of  thermal  stratification  is  to 
be  discussed.  Application  of  the  present  approach  is  being  made  to 
study  the  formation  and  dissipation  of  marine  fog  over  the  ocean. 

2 Analysis 

The  basic  equations,  which  describe  the  atmospheric  boundary 
layer,  consist  of  continuity,  momentum,  temperature,  and  concen 
t ration  For  two-dimensional,  steady  state  boundary  layer  flows, 
the  governing  equations  may  be  written  as 
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and  Fm,  Fh,  and  Fcm  are  the  source  terms  of  momentum,  energy, 
and  concentration,  respectively.  Detailed  discussion  of  the  source 
terms  is  given  in  the  Appendix.  Closure  of  this  system  of  equations 
is  accomplished  by  using  the  Boussinesq  (I8j  approximation  and 
the  Nevzglajdov  (12)  model,  as  modified  by  Lee  and  Harsha  |14j. 
the  apparent  turbulence  kinetic  energy,  Q,  can  be  related  to  the 
laminar  and  the  turbulent  shear  stresses  as 


df 


n * *r  /•  ar 

0.3  - 


It  is  to  he  noted  that  the  magnitude  of  (J  is  contributed  I 
shear  stresses  of  the  eddy  motions  in  the  turbulent  layer,  the 
lecular  motions  in  the  laminar  sublayer,  and  both  laminar  and 
bulent  motions  in  the  transition  layer  An  additional  gover 
equation  for  the  apparent  turbulence  kinetic  energy  is  written 
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concentration  for  species  m.  T is  the  mean  temperature;  the 
primes  denote  fluctuating  quantities,  u.  <*.  and  I)  are  the  dynamic 
viscosity,  thermal  diffusivitv  and  species  diffusivity,  respectively; 


where  Fu  is  the  source  term  for  the  production  and  dissipation  <*t 
turbulence  kinetic  energy. 

Introducing  the  apparent  Prandtl  and  Schmidt  numbers  as  dis- 
cussed by  Lee,  et  al.  (19|.  the  exchange  coefficients  of  momentum, 
heat,  concentration  and  turbulence  become 


A' 


—J,11' 
1 u /~7p 


0.3  pQ  — 

(7) 

4-  Am  'oh 

7 

181 

-Nomenclature* 


n = empirical  constant 
Cm  » concentration  of  species  m 
c'  = concentration  fluctuation 
< > = specific  heat  at  constant  pressure 
(\  = rate  of  formation  of  liquid  water 
f)  = coefficient  of  molecular  diffusion 
/>,,  = dissipation  of  turbulence  kinetic- 
energy 

/•'  "i  = source  term  for  concentration 
equat  ion 

= .source  term  for  temperature 
equation 

F„,  = source  term  for  momentum  ecjua 
tion 

h a source  term  for  turbulence  kinetic- 
energy 

g = giavitat mnul  acceleration 
K = exchange  coefficient  of  concert ra 
tion 
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exchange  coeffic  lent  of  heat 
exchange  coefficient  of  momentum 
exchange  coefficient  of  turbulence 
kinetic  energy 
mass  absorption  coefficient 
latent  heat  of  vaporization 
static  pressure 
pressure  fluctuation 
production  of  turbulence  kinetic 
energy 

kinetic-  energy  of  turbulence 

radiat  ion  energy- 

gas  constant  of  air 

gradient  Richardson  number 

flux  Richardson  number 

mean  temperature* 

temperature  fluc  tuation 

mean  velocity  in  the  .V  direction 

velocity  fluctuation  in  the  X 
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direct  ion 

sedimentation  vt lex  il\  of  fog  drop 
lets 

mean  velocity  in  the  7.  iirection 
velocity  fluctuation  in  the  / direc 
tion 

horizontal  distance 
vertical  distance 
t herrnal  diffuaivit v 
ennssiv  it  v 

boundary  layer  thu  kness 
dv namic  viscosity 
density 

Stefan  Boltzmann  i nst.mt 
turbulent  St  hmnlt  nunibc  r 
turbulent  Prandtl  number 
equivalent  Prandtl  number 
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res|>ectivclv . in  which  «»*  is  the  apparent  Brandtl  numlier  with  a 
numerical  value  of  0.7f>,  obtained  experimentally  by  Arya  (20)  for  a 
simulated  stable  atmosphere,  <r , the  apparent  Schmidt  niiml>er. 
equal  to  ii.7.r>  as  given  by  Fleagle  and  Businger  ( *2 1 ) . and  the 
equivalent  Brandtl  number,  which  is  approximately  equal  to  0.70 
as  determined  by  Lee.  et  al.  | *2*2 1 . Both  Brandt  I and  Schmidt  num 
l>ers  were  allowed  to  vary  from  0.4  to  1.0  under  various  conditirns 
of  thermal  stability  but  proved  to  be  insignificant  in  altering  the 
results  to  any  appreciable  degree.  For  lack  of  any  reliable  function 
al  relation  to  account  for  the  variation  in  the  apparent  Brandtl  and 
Schmidt  numbers,  <ih  and  a,  were  assumed  to  be  constant  in  the 
study 

Using  the  concept  of  exchange  coefficients,  the  governing  equa 
tions  of  momentum,  heat.  mass,  and  turbulence  are  parabolic  type 
The  system  of  equations  can  be  simultaneously  solved  tor  the  de 
pendent  variables  of  velocity,  temperature,  concentration,  and  tur 
bolence  using  the  numerical  procedure  of  Batankar  and  Spalding 
{ 10|  It  should  he  noted  that  the  exchange  coefficients  are  deter 
mined  locally  from  the  apparent  turbulence  kinetic  energy  which  is 
a function  of  the  space  coordinates.  The  detailed  numerical 
scheme  and  all  necessary  boundary  conditions  are  given  by  Bepper 
(1|.  For  the  purpose  of  testing  the  predictability,  the  developed 
method  is  verified  with  wind  tunnel  data  before  being  applied  to 
atmospheric  problems. 

3 Verification  With  Wind  Tunnel  Data 

Verification  of  the  present  model  is  made  by  using  available  ex- 
perimental data  Comparison  with  wind  tunnel  data  is  made  under 
conditions  similar  to  neutral,  stable,  and  unstable  atmospheric 
boundary  layers  While  a number  of  experimental  investigations 
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Fig.  3 Comparison  of  results  for  a stable  atmosphere 


have  been  made  under  neutral  atmospheric  conditions  in  meteoro 
logical  wind  tunnels,  only  a few  cases  contain  turbulence  kinetic 
energy  data  which  are  pertinent  to  the  exchange  coefficient  model 
employed  here.  Similar  investigations  of  experimental  work  in 
thermally  stratified  boundary  layers  disclose  even  less  available 
data 

Breliminary  investigations  are  made  in  a neutral  atmosphere  by 
comparing  the  predicted  results  with  the  wind  tunnel  simulation  of 
Sehon  and  Mery  |lb|  as  shown  in  Fig.  2.  The  subscript,  o,  denotes 
initial  conditions.  A nondimensional  apparent  turbulence  kinetic 
energy  profile,  based  upon  measurements  by  Klebanoff  |l.r>|.  is 
used  to  generate  the  initial  profile  in  the  numerical  scheme  The 
discrepancy  between  measured  and  predicted  turbulence  kinetic 
energy  is  due  to  the  use  of  the  Klebanoff  profile  which  is  for  fully 
developed  flow,  while  the  initial  profile  for  the  Sehon  and  Mery 
case  is  still  a developing  flow  However,  the  necessity  of  predicting 
various  atmospheric  conditions,  without  having  measured  turbu 
lence  data  as  the  initial  profile,  limits  our  choice  to  Klebanoffs 
data  for  all  studied  cases.  The  gins!  agreement  in  velocity  profile 
comparison  indicates  that  a fully  developed  turbulent  profile  for 
overage  velocity  does  not  necessarily  mean  that  the  turbulence  is 
in  equilibrium. 

Investigations  of  stably  stratified  atmosphere  are  made  by  com 
paring  the  predicted  temperature  and  turbulence  profiles  with  the 
wind  tunnel  simulation  data  of  Arva  |20|,  as  shown  in  h ig  < I he 
initial  temperature  distribution  s assumed  to  be  uniform  I he 
Klebanoff  turbulence  kinetic  energy  distribution,  similar  to  that 
for  the  neutral  case,  is  used  for  the  initial  profile  The  underpred 
iction  of  turbulence  kinetic  energy  occurs  within  10  |>ercent  of  the 
boundary  layer  near  the  surface  This  discrepancy  suggests  that 
the  Nevzglajdov  model  is  questionable  in  the  (lose  vicinity  of  the 
wall  boundary,  if  the  measurements  are  truly  reliable  Once  out 
side  this  region  of  the  calculated  maximum  turbulence  kinetu  en 
ergy.  the  predicted  profiles  begin  to  agree  with  the  experimental 
data  Because  the  gradient  Richardson  number,  H,  serves  as  a 
quantitative  measure  of  the  thermal  stratification,  the  distribution 
of  H,  is  plotted  as  a function  of  2/A  in  Fig  4 Although  a discrep 
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Fig.  4 Richardson  number  comparison  for  stable  atmosphere 


am  v occurs  in  the  turbulence  kinetic  energy  comparison,  the  Rich 
anlson  number  predictions  agree  well  with  Arva’s  data  even  in  the 
lower  Ml  percent  of  tin-  boundary  layer.  The  gradient  Richardson 
number.  K is  related  to  the  flux  Richardson  num her,  Jt,.  bv  the 
ratio  of  the  momentum  exchange  coefficient  to  the  heat  exchange 
coefficient. 

Investigations  of  unstably  stratified  atmospheres  are  made  by 
comparing  the  predicted  temperature  and  concentration  profiles 
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Fig  5 Comparison  of  results  for  an  unstable  atmosphara 


Fig.  6 Schematics  of  advection  fog  over  the  ocean 


with  the  wind  tunnel  simulation  data  of  Malhotra  and  Cermak  )2dJ 
as  shown  in  Pig.  5.  Using  the  measurer!  profiles  at  a distance  of  0 d 
m downstream  of  a simulated  source  (Xs  = O.d  m ) as  the  initial 
concentration  and  temperature  profiles,  the  predicted  results 
agree  well  with  the  measured  data  throughout  the  flow  field  No 
turbulence  data  are  available  for  comparison  even  though  the  stan 
dard  Klebanoff  profile  i used  as  the  initial  turbulence  in  the  i.u 
merical  scheme.  Predictions,  with  reasonable  accuracy  tor  turn 
averaged  quantities  of  velocity,  temperature,  and  concentration, 
appear  to  be  possible  using  an  approximate  initial  profile  tor  tur 
bulence;  if  the  turbulence  energy  is  being  balanced  simultaneously 
with  the  momentum,  heat,  and  mass,  consistently  throughout  a 
thermally  stratified  atmospheric  boundary  layer.  Extension  «>l  tins 
method  is  being  made  to  predict  the  formation  of  advection  fog  in- 
curring over  cold  sea  surfaces 

i Advection  Fog  Formation 

Because  of  insufficient  information  regarding  ocean  log-  few 
studies  are  available  in  the  literature  How*  ver.  studies  >n  Ian  1 log 
and  coastal  fog  have  been  made  by  Mack,  et  al  | J li  Fisht  r and 
( apian  (2.r>{.  Zdunkowski  and  Trask  (2h'|,  and  Barker  ■.'■l  < > i>id 
eration  of  the  microphysical  processes  of  condensation  radiation 
and  sedimentation  in  developing  the  advection  log  model  is  based 
on  the  current  information  used  for  land  and  coastal  fog  st, id.es 

K;g.  6 shows  schematically  the  formation  of  advection  fog  ,i- 
warm  moist  air  blows  over  an  aqueous  surface  with  > out u..  . . -U 
decreasing  isotherms.  As  a result  of  the  transfer  of  thermal  energy 
between  air  and  ocean,  the  air  temperature  near  the  ocean  surf;  i * 
is  lowered  until  the  dew  point  temperature  is  reached  Tins  pm* cs* 
results  in  the  condensation  of  water  vapor  into  small  liquid  drof 
lets,  f it  is  assumed  that  suffii  lent  number  of  sea  salt  partu  les  ,ir« 
available  ns  condensation  nuclei  A brief  outline  t th*  satin  a n 
adjustment  procedure,  developed  by  M<  Donald  1 K «n*t  A*».  i n 
is  given  in  the  Appendix 

Fig  7 shows  the  effect  of  wind  velocin  >n  fog  ivu-h  ; u»  .t 
Wind  velocities  of  d and  fi  m s are  used  for  ..  uib  ig  an  ..  m>* 
spheric  boundary  layer  with  a thickness  of  'U  • I i1  ! • 1 

stant  lemperature.  seen  to  he  emitting  from  the  • I <.  h • 

t he  continuous  cooling  of  t he  surfac  e.  as  show  u I or  < n*  w m< 1 • 
tv  of  h m/s  case  A similar  trend  is  also  found  in  iln  wind  i 
of  d m/s  case  Defining  the  fog  boundary  as  iht  limn  where  th*  hq 
uid  water  content  approaches  /.» ro,  it  is  m ted  that  an  in*  >*a  ■ o 
wind  velocity  results  in  a decrease  of  th  t g la\*r  I ■ , If  M*  .• 
over.  fog  only  forms  m the  region  where  a stable  tin  r m.tl  >i  i ,i  du 
turn  occurs  This  result  ran  he  qualitatively  verifw  d front  fn-i  ‘ >h 
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servati«»ns  Leipper  j 30)  observed  that  !*>r  a period  of  two  days,  fog 
formed  very  rapidly  when  the  temperature  distribution  ot  a stable 
atmosphere  incurred  near  the  ocean's  surface.  Liquid  water  con- 
densed out  continuously  to  form  fog  particles  us  long  as  the  stable 
stratification  persisted  Until  the  temperature  distribution  near 
the  ocean's  surface  returned  to  a neutral  atmosphere  on  the  third 
day.  fog  began  to  dissipate  Visibility  started  to  improve  on  the 
fourth  day  as  the  temperature  inversion  region  lifted  up  from  the 
ocean  surface  to  the  lower  atmosphere  Fog  became  stratus  cloud 
ri  the  tilth  day  when  the  lower  atmosphere  returned  to  neutral 
condition. 

5 Conclusion 

An  analytical  investigation  of  turbulent  diffusion  in  thermally 
stratified  boundary  layers  is  made  by  using  a phenomenological 
theory  based  upon  the  apparent  turbulence  kinetic  energy  The  ex 
change  coefficient  of  momentum  is  determined  based  on  the  modi 
tied  Nevzglajdox  model  which  relates  the  local  shear  stress  with 
the  local  turbulence  energy  The  exchange  coefficient  of  heat  and 
mass  are  determined  bv  relating  them  with  the  exchange  coe' 
cient  of  momentum  through  Prandtl  and  Schmidt  numbers.  These 
coefficients  can  only  be  evaluated  simultaneously  with  other  de 
pendent  variables  of  velocity,  temperature,  and  concentration. 
Verification  of  this  method  is  made  by  using  wind  tunnel  simula 
lions  of  thermallv  stratified  atmosphere  Predictions  for  actual 
condition-  can  only  be  quantitatively  verified  when  more  com 
pleted  lielc.  data  are  made  available.  This  study  demonstrates  that 
the  development  of  a realistic  model  for  predicting  the  transport 
phenomena  between  atmospheric  and  aqueous  surface  layers 
needs  close  cooperation  among  researchers  in  many  disciplinary 
areas. 
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Tin  source  terms  used  in  the  governing  equation*  are  derived 
from  inlormntion  avadalile  in  the  current  literature  F<<r  the  read 
i f i oitvetiiem  e.  a hriel  summary  is  given  in  the  following 

Momentum  Equation.  Pressure,  gravity,  buoyant  v and  ( < n 
• lis  f.>r«  i s mav  all  be  involved  in  the  sour»  e term,  however,  in  the 
-t 1 idu-d  . uses,  only  the  pressure  lorce  is  significant 


Fig  7 Fo^iayar  height  affected  by  wind  velocity 


Transactions  ot  the  ASME 


FEBRUARY  1975 


I- m (»P/->X)  (A  - 1 ) 

and 

r P^expt-  X/JllJ)  (A-2) 


in  which  P,i  is  the  atmospheric  pressure  at  sea  level,  and  li„  is  the 
gas  constant  of  air 

Temperature  Kquation.  Latent  heat  and  radiation  are  con 
sidered  by  the  relation 

F,  [ / »(' , -{9R/aZ)\/c,  (A-3) 

in  which  /, n is  the  latent  heat  of  evaporation  at  a given  tempera 
ture.  and  is  the  rate  of  formation  of  condensed  liquid  water  at  a 
given  temperature  and  its  corresponding  partial  pressure  of  vapor 
In  order  to  calculate  the  saturation  adjustment  procedures  of 
McDonald  (28|  and  Asai  (29j  may  he  outlined  as  follows.  Assuming 
air  and  water  vapor  can  he  treated  as  ideal  gases,  the  amount  of 
condensed  liquid  water  can  he  expressed  as 

0.622  P'  ~ 0.622  (A -4) 


for  each  unit  of  air  mass,  in  which  l\,M  is  the  partial  pressure  of 
vapor  at  saturation,  Pv  is  the  partial  pressure  of  vapor  at  each 
thermodynamic  state  of  the  atmosphere,  and  P„  is  the  partial 
pressure  of  air  which  is  approximately  equal  to  the  total  pressure 
of  the  atmosphere.  P By  using  the  Clausius-Clapeyron  equation 
for  liquid  vapor  phase  change,  the  change  of  vapor  pressure  corre- 
sponding to  the  change  of  temperature  can  he  written  as 


T,  - 7 K.,7- 


IA-5) 


The  saturation  vapor  pressure.  PK,  which  is  a function  of  tempera 
ture  only,  is  given  by  Murray  [31 J as 


,17.269(7-  273.16), 

6.108  exp{ T~~~35Ji6 1 (A'6) 


However,  the  first  law  of  thermodynamics  requires  that  :he  pres- 
sure and  temperature  of  vapor  in  an  isobaric  atmosphere  satisfy 
the  following  relation 


l\.  - i\.  <V’ 

T,  - 7 0.6221.,, 


(A -7) 


in  which  rp  is  the  specific  heat  of  air  By  eliminating  l\.  between 
equations  (A-5)  and  (A-7),  the  change  of  temperature  before  and 
after  condensation  becomes 


0.622 (/>  - 

h pKJ-  + 0.387/.^  ' 


(A-8) 


Bv  combining  equation  (A-8).  (A  6).  and  lA-4).  the  amount  of  liq 
uid  water  condensed  per  unit  air  mass  becomes 


J*<T.  ~ T). 


(A -9) 


The  rate  of  formation  of  liquid  water  may  then  he  written  .is 


C,  f pVudZ  (A -10) 

A 

The  radiation  term  occurring  in  equation  (A  3)  may  he  written  in 
the  form  suggested  by  Mack,  et  al.  |23| 


H(Z)  /io7 ' 1 exp  (- 1.6 k\J  pnsmdZ)  A-ll) 

' *b 

in  which  fj  is  the  emissivity  of  a nonhlack  surface  <d  * 0 25),  o is 
the  Stefan  Boltzmann  constant.  7’u  is  the  surface  temperature,  k, 
is  a single  spectrally  averaged  mass  absorption  coefficient  ik,  * 
1.5  X JO  crn/grn).  and  Z,  and  Z/,  are  the  top  and  bottom  of  the  fog 
layer,  respectively  The  gradient  of  the  radiation  effect  then  In- 
comes 

1.6(107  *kj> U.”  exp(-  1.6fcx  f ' ps’ilZ)  (A  - 1 2) 

<)Z  “ 

in  which  u)"'  is  the  liquid  water  content  per  unit  air  mass 

Species  Kquation.  Three  species  are  considered  in  the  »•> 
model:  air  (m  = 1 ).  water  vapor  (m  = 2),  and  liquid  water  i >n 
No  air  can  be  formed  during  fog  formation  and  dissip/i in  • pr 
cesses.  Consequently, 

F,.'  0 'AH 


Water  vapor  can  only  be  formed  at  the  expense  t l - • , • ' 
thus 

Fc  -C\  iA- 1 4 

Liquid  water  can  be  formed  by  condensation  of  water  vaj »>r  and 
also  be  removed  by  fallout  as  water  drops  This  leads  to  the  expre- 
sion 

Fe%  < , -*  ~ (A -15) 

where  Vd  is  the  falling  velocity  of  water  drops 

l'd  5.3  • 10  )“  (A -16) 

with  N representing  the  numiter  ol  drops  par  unit  volume  Jt  is 
noted  that  equation  (A  18)  is  based  on  the  assumption  that  all 
drops  are  less  than  20  n\r  in  radius.  Assuming  a drop  concentre 
tion  of  N — 50  cm*1, 

rd  4 * 102(u>m)2n  (cm/s).  (A -17) 

Turbulence  Kinetic  Energy.  Production,  Pw,  and  dissipation. 
I)u.  of  the  turbulence  kinetic  energy  are  discussed  by  Lee.  et  al 
|22|  and  Plate  |32|.  'The  production  term  may  be  written  as 


which  consists  of  turbulence  energy  generated  by  the  mean  shear 
How  and  the  thermal  stratification.  The  dissipation  term,  as  dis 
cussed  by  Patankar  and  Spalding  (10)  and  Clushko  |9).  may  be 
written  as 

l)„  (A -19) 


in  which  0 is  the  boundary  layer  thickness.  The  coefficient  for  non 
isotropy  of  the  turbulence,  a.  is  given  bv  Byrne  and  Lee  |33[  as 


(i  1.8  for  / 

A 

" u,r/~  4 


(A  -20) 


whirh  agrees  with  the  model  used  hy  hradahaw.  el  al  (H| 
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Numerical  Methods  for  Separated  Flow  Solutions 
around  a Circular  Cylinder 


C.  I . l in,*  I).  W.  Pepper, t and  S.  C.  Lec+ 
University  of  Missouri- Rolla,  Kolia,  Mo. 


Numerical  solutions  of  the  Navier-Slokes  equations  were  obtained  tor  separated  Hows  around  a circular  cylin- 
der at  Reynolds  numbers  40.  HO.  and  2(H).  I he  flow  fields  were  obtained  b>  using  three  finite-difference 
techniques.  I he  implicit  scheme  solved  h>  matrix  factorizations  gave  the  best  accuracy  and  used  the  least  com- 
puter time  I he  flow  pattern  hi  the  recirculating  region  of  a circular  cylinder  begins  to  oscillate  as  the  Reynolds 
number  exceeds  40.  I he  calculated  drag  coefficients,  separation  angles,  and  Slrouhal  numbers  were  compared 
with  available  experimental  data.  < omputational  inaccuracy  resulting  from  numerical  approximations  needs  to 
he  identified  before  a complicated  flow  phenomenon  n he  realistically  analyzed. 


Nomenclature 

a arbitrary  constant  tor  coordinate  transformation 

/ - frequency  of  Karman  vortex  street 

l/.J  lower  matrix  after  modification 

\\f]  = coefficient  matrix 

rti  — number  of  iterations 

|.V]  modifier  of  the  coefficient  matrix 

n = number  of  time-steps 

\q\  - column  matrix 

r - radial  coordinate 

K radius  of  the  cylinder 

He  Reynolds  number,  2UR/v 

S - Slrouhal  number,  2/R/U 

t =time 

[(.]  upper  matrix  after  modification 

U = freest  ream  velocity 

l'„  = velocity  in  the  0 direction 

\ , = velocity  in  the  r direction 

/ = nondimensionali/ed  radial  coordinate 

w ) =vorticity 

\l>  = stream  function 

i - kinematic  viscosity 

0 = angular  coordinate 

At  = time  increment 

At/  = angular  increment 

A/  radial  increment 

1 4*  | :olumn  matrix 

I.  Introduction 

BLC  \l  SI  the  Nax ter  Stokes  equations  are  non- 
linear, exact  solutions  arc  not  currently  available.  The 
necessity  of  providing  reasonable  estimates  for  complicated 
flow  phenomena  leaves  research  engineers  xery  little  choice. 
The  numerical  method  is  one  of  the  very  few  acceptable  tools 
that  is  capable  of  making  any  contribution  to  engineering 
designs  or  environmental  controls  Recent  developments  in 
high  speed  digital  computers  make  this  approach  both  cl 
fectivc  and  popular,  however,  every  numerical  approach  for  a 
given  physical  problem  requires  both  physical  assumptions 
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and  mathematical  approximations  Computer  solutions, 
especially  in  practical  problems  of  questionable  physical 
assumptions  such  as  an  unverified  closure  scheme  for  tur 
buleni  flows  as  discussed  by  1 ee  and  Harsha  and  Pepper  and 
Lee.:  may  appear  to  give  reasonable  results  by  substituting 
unrealistic  physical  assumptions  with  incorrect  mathematical 
approximations.  One  needs  to  test  a numerical  method  for  its 
accuracy  in  the  developing  stage  by  applying  it  to  a given 
problem  that  requires  no  physical  assumptions.  Separated 
flow  around  a circular  cylinder  at  relatively  low  Reynolds 
numbers  is  one  of  the  thoroughly  investigated  problems  that 
can  provide  the  basis  for  such  a study  in  numerical  accuracy 
Numerical  techniques  for  solving  partial  differential 
equations  have  been  discussed  by  Conte.  Richtnner. J 
Roache.  and  many  others.  Solutions  for  separated  flow 
problems  usually  require  substantial  computer  time  because 
of  the  elliptical  nature  of  the  governing  equations 
Multidimensional  separated  flows  occur  in  many  engineering 
and  environmental  problems,  which  require  not  only  the 
solution  of  the  Nax ier  Stokes  equations  but  also  the 
simultaneous  solutions  of  the  energy  and  species  equations 
Past  experience  indicates  that  some  numerical  techniques  may 
give  better  accuracy,  whereas  others  may  need  less  computer 
time.  To  avoid  an  unnecessary  waste  of  effort  and  resources, 
we  conducted  a preliminary  investigation  in  which  we  applied 
several  commonly  used  numerical  methods  to  a thoroughly  in 
vestigated  problem  for  the  purpose  of  comparing  the 
numerical  accuracy  and  the  required  computer  time. 

I he  problem  of  separated  flow  around  a circular  cylindet 
has  been  thoroughly  investigated  both  experimentally  and 
theoretically.  Available  information  on  drag  coefficients, 
separation  angles,  and  Slrouhal  numbers  has  been  reported  by 
Schlichting/'  Thoman  and  Szewc/yk,  and  many  others  In 
the  present  study,  physical  evidence  will  be  used  lor  m 
vestigating  the  effect  of  numerical  approximations  on 
solutions  of  elliptical  partial  differential  equations. 

II.  Analysis 

The  governing  equations  for  an  incompressible,  two 
dimensional,  unsteady  laminar  flow  over  a circular  cylinder 
can  be  written  as 

do)  l a ( zw  \ . ) / I I / d I do:  \ 
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■ :>  • ^ • ' ; "/  >'r'’  • 


1 r 


:¥?B3 


Jll  Y 1 97ft 


SEPARATED  FLOW  SOLUTIONS  AROUND  A CIRCULAR  CYLINDER 


SOI 


In  these  equalions,  ui  is  the  sort  icily , and  I and  L.,  are  the 
velocity  components  in  the  rand  ft  directions,  respectively.  In 
terms  of  the  stream  function  4 the  velocity  components  may 
be  sv  r it  ten  as 


V,  = ((l/r)(94/90)\.  V.,  = 94/9r 


(3) 


When  radius  R is  used  as  the  characteristic  length  and  the 
treestream  velocity  U as  the  characteristic  velocity,  the  non- 
dimcnsionali/ed  stream  function  and  vorticity  become 


4 — 4 / (2R , uJ  — u .'R  / U 

The  independent  variables  of  r,  r,  and  flare 

f =tU/R.  r'  = r/R,  fl’=8/a 


(4) 


(5) 


in  which  a is  an  arbitrary  constant  for  controlling  the  in- 
crement ol  the  transformed  coordinate.  Because  the  flowfield 
variations  take  place  more  rapidly  in  the  vicinity  of  the  cylin 
der  than  in  regions  at  large  distances  from  the  cylinder,  it  is 
convenient  to  transform  the  radial  coordinate  by  an  ex- 
ponential function. 


ej/  = r/R 


(6) 


I he  nondimensionali/ed  governing  equalions  with  the  primes 
eliminated  for  simplicity  then  become: 
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ui  which  Re  is  the  Reynolds  number  based  on  the  diameter  of 
the  circular  cylinder  and  the  treestream  velocity. 

III.  Finite-Difference  Methods 

Numerical  solutions  were  obtained  by  writing  the  governing 
differential  equations  of  vorticity,  Eq.  (7),  and  stream  func- 
tion, Eq.  (K),  into  finite  difference  forms.  By  using  i and  j to 
denote  the  locations  in  the  fl  and  7 directions,  respectively, 
and  employing  m for  the  number  of  iterations  and  n for  the 
number  of  lime-steps,  numerical  solutions  for  the  separated 
flow  around  a circular  cylinder  were  obtained  with  the 
following  three  methods. 

\ lillf-t.sl 

flic  directional  difference  explicit  (DDE)  method,  as 
discussed  by  Thoman  and  S/ewc/yk,  was  used  for  solving 
the  vorticity  equation.  The  values  of  w, , at  the  ( n -t-  / ) th  time- 
step  were  calculated  hv  using  those  of  the  (rrlth  lime  step 
through  the  following  relation 
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Once  the  vorticity  values  were  obtained,  the  stream  functions 
could  be  calculated  by  using  the  Gauss-Seidel  iteration  (GSli 
method  as  discussed  by  Conte. J The  values  of  4 ai  the 
(n  + 7)(h  time-step  and  the  (m  + / ) th  iteration  were  calculated 
by  using  the  following  relation 

a ”i  • i - 

<7 
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(11) 


in  which  d, j = [2/  (A/)  * j + \2 / (AO)  T he  iteration  process 
was  completed  vs  hen  the  difference  of  the  numerical  values  of 
\l  between  the  (m  + /)th  and  the  (w)th  iterations  vsas  less 
than  a prescribed  tolerance.  The  same  procedure  was  repeated 
for  the  successive  time-steps  until  steady-state  solutions  were 
reached . 

B.  ADI-SOK 

The  alternating  directional  implicit  (ADI)  method  was 
proposed  by  Peaceman  and  Rachford * for  solving  linear  dif 
ferential  equations  and  modified  bv  Douglas,1*  Wachspress. 
and  many  others  for  improving  computational  accuracy  Son 
and  Manratty"  used  the  method  for  flow  around  circular 
cylinders  with  the  assumption  that  the  How  field  is  symmetric 
to  the  axis  of  the  How  direction.  I in  and  Lee1'  used  the 
method  for  calculating  transient  flows  around  a sphere 
However,  it  was  noted  that  for  most  Reynolds  numbers 
separated  flow  oscillates  in  the  wake  region.  Bv  eliminating 
the  assumption  of  flow  symmetrv  in  the  wake  region,  an  ex 
tension  of  I in  and  I ec\  approach  was  used  for  solving  the 
vorticity  equation.  I wo  half-time  steps  were  used  The  lust 
hall  i me  tep  is  an  implicit  finite  difference  equation  in  the  0 
direction 
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I quatmn  Cl)  is  an  explicit  form  of  the  nonliticat  vorticitv 
equation  In  order  to  maintain  numerical  stability,  the 
follow  mg  icslrietions  were  suggested  bv  I human  and 
S/ewc/vk  lor  the  nonlinear  terms 
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I he  second  half  lime  step  is  an  implicit  finite-difference 
equation  m the  / direction 
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I he  SIP  method  uses  a technique  of  matrix  factorization  and 
elimination,  as  discussed  by  Stone14  and  Weinstein  et 
al  r Detailed  procedures  are  given  in  the  Appendix  The 
values  of  stream  function  were  obtained  by  the  same  SIP 
method  through  the  application  of  the  following  finite 
difference  equation 
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Solutions  of  u>,  at  the  (n  + 1 )th  time-step  were  used  in  a suc- 
cessive over -relaxation  (SOR)  method  for  calculating  the 
stream  functions  by  the  following  relation 

<1 = i£"  + * [a -S',j  i.itD! 

■ -.V.rytfr./J  db) 

in  which 

/ j.  = 2 2 

(AZ)’’  (AO)’’  (A/)-'  (AH)’’ 

s . , , =S,  S;  . , =.S,;  , (I7a-e) 

The  term  A indicates  the  optimum  relaxation  factor.  The  value 
ol  A,  which  equals  1 .65,  was  used  in  the  analysis  based  on  the 
method  given  by  Carre.'1  The  stream  functions  \l  at 
the  in  + / )t h time  step  were  obtained  when  the  difference  be- 
tween the  (m  + /)th  iteration  and  the  (m)th  iteration  was 
within  a prescribed  tolerance.  The  same  procedures  were 
repeated  for  the  successive  time-steps  until  the  steady  state 
was  reached. 


r;  =D  , H'  - H'  . q =aV,,/w\7  7 (21  a-f) 

Theoretically,  the  SIP-SIP  method  could  be  used  directly  for 
steady-state  solutions  if  the  boundary  conditions  were  known 
Our  problem  was  that  the  vorticity  values  at  the  surface  of  the 
cylinder  could  not  be  prescribed.  Moreover,  we  were  in- 
terested in  investigating  the  oscillation  phenomenon  in  the 
wake  region.  Consequently,  vorticity  solutions  of  successive 
time  s.eps  were  necessary. 

IV.  Boundary  Conditions 

Based  on  no-slip  conditions,  the  velocity  components  ate 
zero  at  the  surface  ol  the  cylinder,  whereas  a uniform 
How  field  of  velocity  U surrounds  the  cylinder.  The  boundary 
conditions  for  the  stream  functions  and  vorticities  in  terms  ot 
the  transformed  coordinates  (Z,0)  are  given  in  I ig  1.  It  is 
noted  that  the  vorticity  values  at  the  cylinder  surface  can  only 
be  determined  from  the  second  derivatives  of  the  local  stream 
functions  by  the  relation 

I / 3\  1 

uj=  , for  all  3 (22) 

I a - 3/  I / 

In  finite-difference  form  at  7 = 0,  the  vorticity  values  for  all  H 
become 


( . Ml*  SIP 

The  strongly  implicit  procedure  (SIP)  was  proposed  by 
Stone  4 for  solving  linear  elliptical  partial  differential 
equations.  Solutions  of  vorticity  and  stream  function  were  oh 
tamed  separately  at  each  time-step.  The  vorticity  values  at  the 
(nil)  time  step  were  calculated  by  using  the  following  finite 
difference  equation 

fi  *■  l ) V*  - -f-  / u.’  + / , u v , / -f  / / , u,’ , , , / ~ q 

(18) 

in  w hich 


= W,., . / — 2^,,  *•  . I ' a (A/) : (23) 

On  the  surface  of  the  cylinder,  / - 0,  the  stream  I unction  is 
zero,  1^  , -0.  The  stream  function  inside  the  cylinder. 

, /,  is  assumed  to  be  r he  mirror  image  of  the  stream  fun c 

tion  around  the  cylinder,  v • . because  of  zero  gradient 
Consequently,  the  values  of  ^ , on  t Tic  surface  for  all  H can  be 
written  as 

a*  - 2v  . /a ;( A/)  (24) 
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The  values  of  v . • which  arc  not  known  at  steady  state,  can, 
(19a)  in  an  iteration  process,  be  arbitrarily  assumed  to  determine 
the  steads  state  condition.  The  number  of  iterations,  it  thev 
converge,  may  depend  on  the  initial  assumptions  i hat  require 
(I9h)  some  criteria  in  order  to  have  any  consistency  Because  wc 
were  interested  in  finding  steady-state  solution  for  both 
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oscillating  anti  nonoscillating  flows,  it  was  to  our  advantage 
to  use  the  impulsively  started  flow  as  the  initial  condition.  The 
number  of  time-steps  (NT)  indicates  the  history  of  flow 
oscillations.  Because  the  oscillation  phenomenon  is  a part  of 
two-dimensional  separated  flows,  steady  state  is  defined  by 
the  following  conditions:  I)  For  nonoscillating  flows,  the 
values  of  4>  and  u>  remain  constant  at  a given  Z and  0 location 
as  NT  increases;  2)  For  oscillating  flows,  the  variations  of  ^ 
and  u>  in  the  Z and  ft  plane  repeat  themselves  with  a definite 
pattern  as  NT  increases.  Numerical  solutions  for  separated 
flow  around  a circular  cylinder  were  calculated  at  three 
Reynolds  numbers  (40,  80,  and  200)  by  using  three  different 
numerical  methods. 


as 

_ °°5~^>e  Q0B 


^ QOS JJQ> 


4^-005 

Jl  005 


°°5  ? ^ 05 

005  , J°  -.05 


Hr.  2 Streamline  patterns.  Re  ^ 40. 


V.  Results  and  Discussion 

The  effect  of  mesh  sizes  and  outer  boundary  locations  on 
numerical  accuracy  has  been  discussed  by  Lin  and  Lee12  for 
solutions  of  Navier-Stokes  equations.  The  current  study  was 
undertaken  to  investigate  the  relative  accuracy  and  computer 
time  of  three  different  methods.  Only  one  standard  mesh  size 
of  A/  0.1,  = 6°,  and  A/  0.02  was  used  for  a flowfield 

with  an  outer  boundary  of  90  radii  of  the  cylinder.  To  deter- 
mine the  relative  accuracy  of  each  method,  calculated  results 
of  drag  coefficient,  separation  angle,  and  Strouhal  number 
were  compared  with  available  experimental  data.  It  was 
found  that  the  SIP-SIP  method  gave  the  same  accuracy  as  the 
ADI  SOR  method,  which  proved  better  accuracy  than  the 
DDL  CiSI  method.  The  required  computer  time  on  an  IBM 
370-168  computer  is  shown  in  Table  1.  It  is  evident  that  the 
SIP-SIP  method  is  a more  efficient  method  for  the  in- 
vestigated problem.  In  order  to  provide  a complete  picture  of 
the  investigated  result,  the  discussion  is  divided  into  two 
parts. 

•\.  How  Patterns 

I he  stream  function  pattern  calculated  from  the  SIP  SIP 
method  is  shown  in  Pig.  2 for  the  case  of  Reynolds  number  40. 
If  is  noted  that  the  separated  flow  in  the  wake  of  a circular 
cylinder  is  symmetric  and  steady  as  long  as  the  Reynolds  num- 
ber is  less  than  40.  The  wake  length  is  two  to  three  times  as 
great  as  the  cylinder  diameter,  and  the  separation  angle  is 
about  126°  from  the  leading  edge.  Figure  3 shows  that 
streamline  pattern  for  the  case  of  Reynolds  number  80.  The 
separation  angle  is  about  114°.  The  recirculating  flow  in  the 
wake  region  oscillates  with  a definite  frequency.  The  flow  pat- 
tern appears  to  repeat  itself  between  2500  NT  and  3150  NT 
By  using  the  dimensionless  time  increment  A/  0.02,  the 

I able  I ( ompuler  lime  comparisons 


IBM  370-168 


He  I one  (mini 


025  — 

'Cq C''cC'ar)  ) / / 


/ Ub 

^ 

_ -Q8 

-~-^Q25 

^ / f 


fig.  3 Streamline  patterns, 3 e 80. 

oscillating  frequency  can  be  calculated  by  the  following 
equation 

f^U/RM  l (NT)..  (NT)  1 (25) 

T he  Strouhal  number,  given  b\  the  relation 

S = 2Rf/U  (26) 

has  a value  of  0.154  for  the  case  of  Reynolds  number  80 
Figure  4 shows  the  streamline  pattern  for  the  case  of  Rev  nolds 
number  200.  The  separation  angle  is  about  102P  from  the 
stagnation  point  of  the  leading  edge.  The  How  pattern  appears 
to  repeat  itself  between  2400  NT  and  2950  NT . T his  gives  a 
Strouhal  number  of  0.182  for  the  case  of  Reynolds  number 
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200.  By  using  the  AD1-SOR  method,  the  flow  pattern 
calculations  were  obtained  for  cases  of  Reynolds  number  40 
and  80  only.  The  appearances  of  streamline  distributions  are 
very  much  similar  to  those  of  the  shown  results  which  were 
obtained  by  the  SIP-SIP  method.  The  streamline  patterns  ob- 
tained by  the  DDE-GSI  method  deviated  noticeably  from  the 
shown  results.  In  order  to  answer  the  question  of  the  relative 
accuracy  of  each  method,  it  was  necessary  to  compare  the 
results  with  available  experimental  data. 

H.  < umparisun  wiili  1 xperimenlal  Data 

Experimental  data  are  available  iri  terms  of  drag  coef- 
ficient, separation  angle,  and  Strouhal  number.  The  total 
drag  coefficient  is  calculated  by  the  summation  of  the 
pressure  drag  coefficient. 
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,co,'-=R,nrsL  *-r. « i 

(27) 

and  the  skin  friction  drag  coefficient, 

y i* i 

((,,)  = ~ \ cl* I sin0d0  (28) 

By  using  the  obtained  values  of  vorticity  and  vorticity 
gradient  along  the  surface  ol  the  cylinder,  the  drag  coefficient 
tor  each  Reynolds  number  case  was  determined  and  compared 
with  available  experimental  data. 

Numerical  solutions  were  also  obtained  by  other  in- 
vestigators who  used  one  of  t he  three  different  methods  The 
DDI  ( iSl  method  was  used  bv  fhoman  and  S/cwc/yk  who 
compared  their  results  with  the  experimental  data  ol  Kell  and 
Simmons  as  well  as  with  the  data  ot  Morkovin.  I he  ADI  SOR 
method  was  used  bs  Son  and  Manratty  who  also  discussed  the 
numerical  results  of  kawaguati  and  Jain  When  we  used  the 
1)1)1  (.SI  method  tor  our  calculations  in  the  drag  coefficient, 
we  obtained  practically  the  same  results  as  Thoman  and 
S/ewc/vk;  however,  when  we  used  the  ADI  SOR  method  for 
our  drag  calculations,  we  got  the  same  results  as  those  ob 
tar ned  by  if  SIP  SIP  method,  which  is  developed  in  the 
present  work  I he  assumption  that  a line  of  symmetry  exists 
iri  the  wake  region  appears  to  he  the  source  of  inaccuracy  in 
the  numerical  results  of  Son  and  Manratty  at  Reynolds  num 
hers  greater  than  4<)  When  compared  with  the  experimental 
data  ot  Perry  ' and  Roshko  as  shown  m f ig  5,  this  ob 


fig.  6 Separation  angle  comparisons. 

servation  was  confirmed.  It  is  interesting  to  note  that  the 
discrepancy  between  experimental  data  is  greater  than  that  of 
the  numerical  results  at  the  lower  Reynolds  number  region 
The  experimental  results  were  also  compared  with  analytical 
results  on  separation  angles  as  shown  in  Fig.  6.  The  ex- 
perimental data  obtained  by  Homannlh  at  Reynolds  number 
40  agree  well  with  all  analytical  results.  The  DDE-GSI  method 
of  Thoman  and  Szewc/yk  appears  to  give  the  most  inaccurate 
predictions  at  other  Reynolds  numbers.  I his  is  to  be  expected, 
because  the  criteria  given  in  Eq.  (10)  lor  numerical  stability 
tend  to  generate  larger  values  of  inaccuracy  in  the  close 
vicinity  of  a solid  boundary  where  the  flow  separation  is 
originated.  Moreover,  the  assumption  of  symmetry  in  the 
wake  region  that  was  used  by  Son  and  Hanrattv  apparently 
gives  reasonable  predictions  for  separation  angles  because 
they  are  located  upstream  of  the  wake  region  I he  ex 
penmental  and  numerical  results  on  Strouhal  numbers  can 
also  be  compared  as  shown  in  fig  7.  Because  the  Strouhal 
number  is  related  to  the  shedding  of  vorncities  in  the  wake 
region,  as  given  in  Eqs.  (25)  and  (26),  numerical  method**  to 
predict  Strouhal  numbers  cannot  assume  the  line  ot  symmetiv 
in  that  region  The  experimental  data  of  Roshko  appear  to 
be  consistent  at  all  tested  Reynolds  numbers.  The  ex 
penmental  data  of  Relf  and  Simmons  as  well  as  the  data  of 
Morkovin  were  used  by  Thoman  and  Szewc/yk  for  com 
part. son  with  their  numerical  results,  which  were  obtained  by 
the  DDE-GSI  method.  The  ADI  SOR  method  used  bv  Son 
and  Manratty  cannot  predict  Strouhal  numbers,  because  it 
assumes  symmetry  in  the  wake  region  I he  ADI  SOR  method 
used  m the  present  study  gives  similar  results  as  the  MB  SIP 
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method  tor  the  case  of  Reynolds  number  80  and  uses  more 
than  twice  the  amount  of  computer  time  than  the  latter 
method  (.  onscquently . only  the  SIP-SIP  method  lias  been 
used  tor  studying  the  case  of  Reynolds  number  200.  It  is 
evident  that  discrepancies  exist  in  both  experimental  and 
numerical  results  for  Strouhal  numbers.  The  present  study  ap- 
pears  to  agree  better  with  Roshko’s  data  than  w ith  the  data  in 
other  studies. 

VI.  Conclusions 

It  is  necessary  to  examine  a newly  developed  numerical 
method  for  the  purpose  of  correcting  questionable 
mathematical  approximations.  For  solving  nonlinear  ellip- 
tical differential  equations,  the  numerical  method  which  ap- 
pears to  give  the  best  accuracy  and  uses  the  least  amount  of 
computer  time  is  the  implicit  finite  difference  scheme  that  is 
solved  by  matrix  factorizations. 

Separated  How  around  a circular  cylinder  can  have  a line  of 
symmetry  in  the  wake  region  only  at  Reynolds  numbers  less 
than  40.  In  the  case  of  oscillating  flows,  there  are  two  factors 
which  cause  errors  in  numerical  solutions.  One  is  the 
mathematical  approximation  of  the  nonlinear  terms  that  is 
used  to  achieve  numerical  viability.  The  other  is  the  physical 
assumption  of  the  wake  symmetry.  Calculations  of  the  drag 
coefficient  appear  to  be  within  the  limit  of  experimental  ac- 
curacy it  the  physical  reality  of  flow  oscillations  is  allowed  in 
the  wake  region.  Calculations  of  the  separation  angle  appear 
to  be  affected  more  significantly  by  the  mathematical  ap 
proximation  than  the  physical  assumption.  Calculations  of 
the  Strouhal  number,  however,  require  an  adequate  treatment 
Bpf  both  factors. 

1 able  2 Hementsof 


The  question  on  computer  simulations  vs  wind  tunnel 
simulations  can  only  be  answered  according  to  the  purpose  of 
applications. 

1)  I or  laminar  Hows  at  relatively  low  Reynolds  numbers,  it 
is  economically  advantageous  to  use  computer  simulations 
because  the  technology  available  today  is  adequate. 

2)  As  the  Reynolds  numbers  increase  to  the  region  ot  in- 
stability. computer  simulations  using  unsteady -state  Navicr- 
Stokes  equations  are  questionable,  because  the  numerical  in- 
stability is  not  distinguishable  from  flow  instability. 
However,  wind-tunnel  simulations  may  have  the  same 
problem  of  super-positioning  mechanical  vibrations  upon 
flow  instabilities. 

3)  In  the  turbulent  flow  region,  computer  simulations  need 
to  use  time-steps  at  least  one  order  of  magnitude  smaller  than 
the  smallest  period  of  the  turbulent  fluctuation.  I his  will 
enable  us  to  solve  the  unsteady  state  Navier-Stokcs  equations 
by  consideiing  the  fluctuating  velocities  as  time  dependent 
variables  of  the  momentary  velocity  components.  However, 
today’s  computer  technology  cannot  reach  the  required  coni 
putational  speed.  Computer  simulations  of  turbulent  flow 
problems  are  present (y  being  conducted  by  using  a closure 
scheme  to  relate  the  turbulent  fluctuations  with  the  time 
averaged  components.  Owing  to  the  limited  availability  ot 
reliable  experimental  data,  there  is  a possibility  of  super 
positioning  questionable  mathematical  approximations  upon 
unrealistic  physical  assumptions.  This  is  particularly  true 
when  the  solution  ol  nonlinear  elliptical  partial  differential 
equations  is  required  for  turbulent  flows  in  the  separation 
region.  The  present  study  is  intended  to  minimize  the 
numerical  inaccuracy  before  a closure  scheme  is  introduced 
into  turbulent  flow  problems.  Consequently,  compiler  and 
wind-tunnel  simulations  are  both  necessary  to  improve  our 
current  understanding  of  turbulent  flows. 

4)  Turbulent  flows  with  heat,  mass,  and  momentum 
transfer  often  simultaneously  occur  in  many  practical 
problems,  such  as  combustion  engineering,  environmental 
pollution,  and  weather  prediction.  A realistic  computer 
simulation  technique  will  require  the  advanced  technology  in 
both  computational  speed  and  computer  storage  which  are  yet 
to  be  developed  to  the  operational  stage,  in  the  foreseeable 
future,  wind-tunnel  simulation  is  considered  to  be  the  inoM 
reliable  method  before  full-scale  testing  for  any  successful 
engineering  design  can  be  realized.  However,  the  eeonnmiL.il 
reality  in  conducting  extensive  experimental  investigations 
makes  it  necessary  for  computer  simulations  ot  approximate 
nature  to  be  used  to  supplement  both  wind-tunnel  simulations 
and  full-scale  testing.  For  practical  purposes,  both  computet s 
and  wind  tunnels  are  necessary  to  supplement  each  other  to 
provide  the  reliablily  and  economy  that  are  required  tor  all 
successful  engineering  endeavors. 
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Appendix 

l Miiu  iht-  N|P  Method  lo  Obtain  the  Matrix  Solution 

The  governing  equations,  in  matrix  form,  with  five  nonzero 
diagonal  elements  as  the  coefficient  matrix  may  be  written  for 
either  Eq.  (1 8) or  Eq,  (20) as  follows 
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Lhe  term  (A/]  is  the  coefficient  matrix  of  the  column  matrix 
IT  I , which  represents  either  the  stream  function  C or  the 
vorticity  u>  at  the  («  + / ) th  time-step  for  every  grid  point  in 
the  flow  field.  The  column  matrix  | q"  \ results  from  the 
know  n values  of  the  flow  field  at  the  (n)th  time-step.  The 
elements  of  matrices  (AT],  I <t>  I , and  | q | are  given  in  Table  2. 
In  order  to  accelerate  the  solution  procedure,  matrix  [M]  is 
modified  as  matrix  | V/  + .V).  which  consists  of  seven  nonzero 
diagonal  elements.  Equation  (Al ) then  becomes 
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in  which  rn  is  the  number  of  iterations  for  determining  the 
column  matrix  J «t>  J at  the  (n  + /)th  lime-step.  The  modifier 
matrix  |,V]  has  to  fulfill  the  requirement  that  the  coefficient 
matrix  (A/  + .V]  can  be  factored  into  the  product  of  a lower 
matrix  (/  | and  an  upper  matrix  ((/);  each  consists  of  three 
nonzero  diagonal  elements  in  the  lower  and  upper  portions, 
respectively,  as  follows: 
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The  iteration  parameter  a is  allowed  to  vary  between  zero 
and  unity.  According  to  Stone  ia  values  near  unity  tend  to 
reduce  the  low-frequency  errors  more  rapidly  In  this  analysis, 
ten  values  of  a between  zero  and  unity  are  evenly  distributed 
in  the  range  between  the  high  and  low  frequencies 
The  purpose  of  this  modification  is  to  determine  the  column 
matrix  by  having  the  difference  between  the 

(m+/)th  and  the  (m)th  iterations  within  a prescribed 
tolerance.  If  a column  matrix  | AT  | at  the  (m  + / ) th  iteration 
is  defined  as 


| A4> " ' ' | ( <i> — 1 4> " * ' I"' 
then  Eq.  ( A2)  can  be  w ritten  as 

[M  + N\  | A*  " * ' | ™ * ' = | q " | - (A/)  | <J>  ” * ' | ' 
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If  one  uses  the  lower  and  upper  matrices  given  by  Eq.  (A3), 
Eq.  (A6)  can  be  expressed  as 
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The  use  of  matrices  (L)  and  [U]  accelerates  the  solution 
procedures  by  successive  eliminations.  Once  the  column 
matri>  1 A4> " ‘ 1 at  the  ( w + / ) th  iteration  is  equal  to  or  less 
than  the  prescribed  tolerance,  the  values  of  {<J>'”/|  at  the 
{rn)\h  iteration  become  the  solutions  of  either  the  vorticity  or 
the  stream  function  at  the  ( n + 1 ) th  time-step. 
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